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Abstract 


A discrete random medium is an object in the form of a finite volume of a vacuum or a homoge- 
neous material medium filled with quasi-randomly and quasi-uniformly distributed discrete mac- 
roscopic impurities called small particles. Such objects are ubiquitous in natural and artificial en- 
vironments. They are often characterized by analyzing theoretically the results of laboratory, in 
situ, or remote-sensing measurements of the scattering of light and other electromagnetic radia- 
tion. Electromagnetic scattering and absorption by particles can also affect the energy budget of a 
discrete random medium and hence various ambient physical and chemical processes. In either 
case electromagnetic scattering must be modeled in terms of appropriate optical observables, 1.e., 
quadratic or bilinear forms in the field that quantify the reading of a relevant optical instrument or 
the electromagnetic energy budget. It is generally believed that time-harmonic Maxwell’s equa- 
tions can accurately describe elastic electromagnetic scattering by macroscopic particulate media 
that change in time much more slowly than the incident electromagnetic field. However, direct 
solutions of these equations for discrete random media had been impracticable until quite recent- 
ly. This has led to a widespread use of various phenomenological approaches in situations when 
their very applicability can be questioned. Recently, however, a new branch of physical optics has 
emerged wherein electromagnetic scattering by discrete and discretely heterogeneous random 
media is modeled directly by using analytical or numerically exact computer solutions of the 
Maxwell equations. Therefore, the main objective of this Report is to formulate the general theo- 
retical framework of electromagnetic scattering by discrete random media rooted in the Maxwell— 
Lorentz electromagnetics and discuss its immediate analytical and numerical consequences. Start- 
ing from the microscopic Maxwell—Lorentz equations, we trace the development of the first- 
principles formalism enabling accurate calculations of monochromatic and quasi-monochromatic 
scattering by static and randomly varying multiparticle groups. We illustrate how this general 
framework can be coupled with state-of-the-art computer solvers of the Maxwell equations and 
applied to direct modeling of electromagnetic scattering by representative random multi-particle 
groups with arbitrary packing densities. This first-principles modeling yields general physical in- 
sights unavailable with phenomenological approaches. We discuss how the first-order-scattering 
approximation, the radiative transfer theory, and the theory of weak localization of electromagnet- 
ic waves can be derived as immediate corollaries of the Maxwell equations for very specific and 
well-defined kinds of particulate medium. These recent developments confirm the mesoscopic 
origin of the radiative transfer, weak localization, and effective-medium regimes and help evalu- 
ate the numerical accuracy of widely used approximate modeling methodologies. 
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1. Introduction 


In this Report we discuss fundamental aspects of the scattering of electromagnetic radiation 
by a discrete random medium (DRM), i.e., an object in the form of a distinct finite volume of a 
vacuum or a homogeneous material medium filled with quasi-randomly and quasi-uniformly dis- 
tributed discrete macroscopic impurities called small particles. The general subject of electromag- 
netic scattering is extremely broad and can hardly be fully covered in a single review, which ne- 
cessitates exercising proper selectivity and a careful delineation of the overall scope of the discus- 
sion. Therefore, the main purpose of the opening section is to introduce, in a somewhat ad hoc 
and qualitative manner, several basic definitions and notions and to explain the main focus of this 
Report. 


1.1. Small particles 


The term “small particles” or just “particles” is ubiquitous in the discipline of light (or, more 
generally, electromagnetic) scattering | 1—35] and even enters the very titles of several specialized 
monographs [1,8,11,13,16,18,20,24—26,28-30,34,35]. However, it may not be straightforward to 
give a universal and unambiguous physical definition of this term. For the purposes of this Report, 
a small particle is defined as a small yet optically macroscopic body. More specifically, a small 
particle is a finite discrete physical body that is “small” (or “tiny” or “minute”’) and yet consists of 
a number of atoms large enough that the body can be characterized by bulk optical constants such 
as the electric permittivity, magnetic permeability, and conductivity. The adjective “discrete” 
means that the body can be thought of as having a distinct macroscopic surface separating it from 
the surrounding host medium and acting as an optical interface between the interior and exterior 
materials with different refractive indices. The distribution of the refractive index inside the parti- 
cle does not need to be homogeneous. 

We will see later that the requirement of being characterized by optical constants appropri- 
ate to bulk matter allows one to define rather unambiguously the minimal permissible size of a 
small particle. It is not as straightforward to define the requisite “smallness” of the particle. More 
often than not, the characterization of being small follows from the human visual perception or 
from the need for an optical or electronic microscope to even see the particle. Often however it is 
more appropriate to refer to the optical size of the particle (or its size parameter), defined as the 
ratio of the circumference of the particle’s smallest circumscribed sphere to the wavelength of the 
incident electromagnetic wave. Then the smallness of a particle may be defined by restricting the 
dimensionless particle size parameter to be a few orders of magnitude or less. For example, 
leaves, birds, or decimeter-sized clumps of ice forming Saturn’s rings do not appear to be small 
particles when looked at by a human eye, but they are small particles from the perspective of 
probing them electromagnetically with a remote decimeter-wavelength radar. The factor that 
makes it convenient and possible to define the smallness of a particle in terms of its size parame- 
ter is the fundamental so-called scale invariance rule of electromagnetic scattering. This rule 
states that all dimensionless scattering and absorption characteristics of a finite object depend only 
on the ratio of the object’s size and the wavelength of radiation but not on their individual values 
[36]. 

Fig. 1 shows that in many cases the above definition of a small particle can be rather une- 
quivocal. Fig. 2 illustrates however that a degree of ambiguity can remain in some cases. Indeed, 
on one hand a fractal soot aggregate with touching components can be considered an individual 
particle when it is suspended in the atmosphere and is widely separated from other atmospheric 
particulates. On the other hand, it can also be considered a composite object consisting of individ- 
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ual particles in the form of soot spherules. The possibility of this and similar ambiguities should 
always be kept in mind. 


1.2. Discrete random media 


Using the above definition of a small particle, a Type-1 discrete random medium (DRM; see 
Fig. 3a) can be defined as a morphologically complex object in the form of an imaginary volume 
V populated by a large number N of small particles in such a way that 


e the spatial distribution of the particles throughout the volume is quasi-random and quasi- 
uniform, and 


e the physical states of the individual particles are statistically independent of each other and of 
the particle positions 


(e.g., [7,12,19]). The physical state of a particle is defined as the combination of the particle’s 
size, morphology (including the spatial distribution of the refractive index), and orientation. A 
defining trait of electromagnetic scattering by a DRM is the absence of speckles in scattering pat- 
terns. 

It is imperative to recognize that at any moment in time, the spatial distribution of particles 
in a multi-particle group is definite rather than random. Therefore, if the group is illuminated by a 
monochromatic or quasi-monochromatic' parallel beam of light then statistical randomness and 
spatial uniformity of a DRM and the requisite absence of speckles can be achieved only over a 
sufficiently long period of time owing to random temporal changes of particle positions. This is 
precisely what happens naturally in a multi-particle group suspended in a gas or a liquid and caus- 
es smooth, speckle-free patterns of electromagnetic scattering. 

In some cases (e.g., in a particulate surface) the particles do not move relative to each other 
and yet, from the perspective of electromagnetic scattering, can often be thought of as forming a 
DRM. This can happen, for example, when particle positions in a group are “maximally random” 
and “maximally uniform’, while the entire multi-particle group is moving relative to the source of 
light and/or the detector. Then even small changes of the source-of-light — multi-particle group 
— detector configuration during the measurement are equivalent to multi-wavelength shifts in 
particle positions and can, in essence, result in a random particulate sample generating speckle- 
free scattering patterns [44,45] (see also Section 1.4 of [25]). Another way to achieve a speckle- 
free regime is to illuminate a fixed quasi-random multi-particle group by incoherent polychro- 
matic and/or uncollimated light. Such scattering scenarios help broaden the notion of a DRM and 
often lead to useful practical applications. 

The volume V hosting the N-particle group can also have a distinct physical boundary S sep- 
arating the finite interior and the infinite exterior space with different refractive indices (Fig. 3b). 
Such a heterogeneous object can be classified as a morphologically complex DRM as well, pro- 
vided that the distribution of particle positions throughout the actual physical volume V is suffi- 
ciently random and uniform. We will refer to such an object as a Type-2 DRM. 

In reality, the spatial distribution of the constituent particles can never be completely ran- 
dom and statistically uniform because the particles are not allowed to overlap and because their 
cumulative volume V,,, (defined as the union of the individual particle volumes) is nonzero. For 


' According to conventional terminology, the qualifier “monochromatic” refers to a purely time- 
harmonic electromagnetic field, while the qualifier “quasi-monochromatic” refers to a time- 
harmonic electromagnetic field subjected to relatively slow quasi-random fluctuations. 
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the same reason the orientations of nonspherical particles cannot be completely independent of 
each other and of the particle coordinates. We will assume however that the particles are distrib- 
uted throughout the (physical or imaginary) volume V as randomly and uniformly as the volume 
packing density p = Vyar i V <1 permits.” In this regard the morphology of a DRM is fundamen- 
tally different from that of fractal-like multi-particle clusters such as those studied, e.g., in [41— 
43,4649] and illustrated in Fig. 2. On a somewhat pedestrian level, the spatial distribution of the 
particles can be considered statistically quasi-uniform if the average number of particles per unit 
volume n(r) is independent of the position vector r over distances of the order of several times 
the average separation between two neighbouring particles. Of course, n(1r) is allowed to change 
over much greater distances. An instructive discussion of various mathematical parameterizations 
of physical disorder can be found in [50]. 


1.3. Why to study electromagnetic scattering by discrete random media? 


Discrete random media are ubiquitous in natural and artificial environments. Typical exam- 
ples are clouds of interstellar dust; the cloud of interplanetary dust in the solar system; dusty at- 
mospheres of comets; particulate planetary rings; clouds in planetary atmospheres; geophysical, 
biomedical, and technical particle suspensions; aerosol particles with numerous inclusions; heter- 
ogeneous polymeric materials; and particulate surfaces (cf. Figs. 4 and 5a—d). Another important 
class of DRMs is represented by technical coatings such as layers of paint [54] (Figs. 5e,f). 

The extreme morphological complexity of the majority of natural and artificial DRMs 
makes their characterization a daunting task. More often than not, one has to infer the micro- and 
macrophysical parameters of a DRM by analyzing theoretically the results of laboratory, in situ, 
or remote-sensing measurements of light and other electromagnetic radiation scattered by the 
medium. Thus the use of electromagnetic scattering as a potent noninvasive characterization 
technique represents a major reason to study this phenomenon. Another major reason has to do 
with the fact that scattering and absorption of electromagnetic radiation by particles can affect 
the energy budget of a volume of DRM and hence various ambient physical and chemical pro- 
cesses. In either case electromagnetic scattering must be described in terms of appropriate optical 
observables, i.e., quadratic or bilinear forms in the electromagnetic field that quantify the reading 
of a relevant optical instrument and/or the electromagnetic energy budget. 


1.4. The general scope of this Report 


The practical solution of optical-characterization and energy-budget problems has the fol- 
lowing four main ingredients: 


e Formulation of appropriate optical observables for a given DRM and a specific type of illumi- 
nation. 


e Theoretical modeling of these observables for a specific DRM (the so-called direct scattering 
problem). 


e Practical measurement of these observables. 


e Solution of the so-called inverse scattering problem, 1.e., finding the physical model of a 


* Note that the volume packing density p of a DRM can vary from essentially zero for a cloud to 
more than 0.5 for a particulate surface. 


DRM that provides the best fit of theoretical simulations of electromagnetic scattering to the 
measurement data. 


To keep the size of this Report manageable, we will focus only on the first two ingredients. 
Specifically, our primary objective is to outline the first-principles theoretical framework of elec- 
tromagnetic scattering by DRMs rooted in the microscopic Maxwell—Lorentz equations and dis- 
cuss its immediate analytical and numerical consequences. 


1.5. The need for first-principles approach 


Until quite recently, theoretical calculations of electromagnetic scattering by a DRM had 
typically been based on ad hoc approaches with poorly known or undefined accuracies and rang- 
es of applicability. Perhaps the most notable examples are the phenomenological’ radiative trans- 
fer theory [7,12,19,55—71] originally developed for sparse turbid media such as clouds, effective- 
medium rules [32,72—78], and the geometric-optics ray-tracing method [79-83]. Even simplistic 
phenomenological approaches such as the Gershun theory of the light field [84] or the Kubelka— 
Munk [85-89] and Hapke [90] theories, have found frequent — even though questionable — use 
(cf. [91,92]). The underlying principles of some of these methodologies can be traced all the way 
back to the classical yet thoroughly outdated work by Pierre Bouguer [93,94], Johann Lambert 
[95], August Beer [96], Eugen von Lommel [97], Orest Khvolson [98], and Arthur Schuster [99] 
(see [100] for an account of the history of the phenomenological radiative transfer theory). The 
basic “physically obvious” premise in many studies of electromagnetic scattering by DRMs has 
been the belief that if the individual far-field scattering properties of each constituent particle are 
known then all scattering properties of the entire DRM can somehow be constructed from those 
of the constituent particles. This generally incorrect assumption is based on the lack of recogni- 
tion that from the perspective of electromagnetics, the entire DRM is a unified scattering target, 
while the only essential consequence of the complex object’s morphology (e.g., being composed 
of what appears to the human eye as discrete units, called particles) is to make the corresponding 
electromagnetic boundary conditions more complicated. 

The main objective of this Report is to expose the fundamental physical nature of the phe- 
nomenon of electromagnetic scattering by a DRM and introduce the general theoretical formal- 
ism enabling first-principles modeling of relevant optical observables. We demonstrate how re- 
cent advances in the development of computer solvers of the macroscopic Maxwell equations 
and the availability of powerful computers and computer clusters have made possible direct 
modeling of electromagnetic scattering by representative random multi-particle groups with arbi- 
trary packing densities. Furthermore, we discuss how the first-order-scattering approximation, 
the radiative transfer theory and the theory of weak localization of electromagnetic waves can be 
derived as immediate corollaries of Maxwell’s electromagnetics for very specific and well- 
defined kinds of DRM. These recent developments have decisively brought the subject of elec- 
tromagnetic energy transport in macroscopic DRMs and their optical characterization into the 
realm of physical optics. In particular, they have helped establish a mesoscopic link between the 
macrophysical regime of radiative transfer, weak localization, and effective-medium approxima- 
tions on one hand and the microscopic Maxwell—Lorenz equations on the other. We make a spe- 


3 A physical theory is called phenomenological if it expresses mathematically the results of ob- 
served phenomena without tracing and clarifying their fundamental origin and significance. Typ- 
ically, the development of a phenomenological theory is based on heuristic (1.e., experience- 
based) shortcuts lacking rigorous justification. 


cial effort to state explicitly what results have been established definitively and what aspects of 
this research discipline necessitate further analysis. 

The unquestionable advantage of the first-principles approach is that it yields the definitive 
physical understanding of the phenomenon of electromagnetic scattering by a DRM and its cor- 
ollaries. However, technical complexities of solving the Maxwell equations directly (both analyt- 
ically and numerically) often diminish the applicability of the first-principles approach to real 
physical systems encountered in practice. As a consequence, it is reasonable to expect that vari- 
ous analytical, phenomenological, and heuristic approximations such as those mentioned above 
will still be widely used in the foreseeable future. Hence an important function of the first- 
principles approach is to characterize the accuracy and range of applicability of these approxima- 
tions. To illustrate this function, we include a discussion of how direct computer solutions of the 
macroscopic Maxwell equations can be used to quantify the errors of such popular modeling 
tools as the first-order-scattering approximation, the radiative transfer equation, the theory of 
weak localization, and the effective-medium approach. 


1.6. Further guidelines 


To make the scope of this Report manageable, we will discuss only elastic scattering of 
electromagnetic waves. In other words, nonlinear optics effects will be excluded by assuming 
that the optical constants of the scattering object as well as of the surrounding medium are inde- 
pendent of the electric and magnetic fields. 

We will also exclude from specific consideration the small Doppler shift of frequency of 
the scattered light relative to that of the incident light due to the movement of particles with re- 
spect to the source of illumination. Furthermore, we will not discuss the scattering of transient 
electromagnetic fields such as ultra-short laser pulses (cf. [101]) and will discuss only frequency- 
domain electromagnetic scattering by assuming that all “quasi-instantaneous’” fields and sources 
are time-harmonic and satisfy the frequency-domain Maxwell equations. In other words, we fo- 
cus on the scattering of a monochromatic or quasi-monochromatic electromagnetic field and as- 
sume that the scattering object varies in time much more slowly than the field. 

In the majority of this Report we will assume that the randomness of a particulate medium 
is ensured by its temporal variability. The extension of the concept of a DRM to a fixed particu- 
late medium illuminated by an incoherent source will be discussed in Section 12. 


2. Electromagnetics, optical observables, and averaging 


The most advanced theory of light—matter interactions available today is quantum electro- 
dynamics (QED) [102-108] followed, in the hierarchy of generality and complexity, by the 
semi-classical approach [109-111] and the Maxwell—Lorentz microscopic electromagnetics 
[112-116]. Since the specific subject of this Report is elastic (i.e., not involving changes in fre- 
quency) electromagnetic scattering, we will assume that from the standpoint of a wide range of 
practical applications, all relevant physics can be adequately captured by the classical microscop- 
ic Maxwell—Lorentz equations. 

Despite this simplification, the actual quantification of electromagnetic scattering by a 
DRM is still highly problematic because solving the Maxwell—Lorentz equations either analyti- 
cally or numerically is essentially impossible given the enormous number of elementary electric 
charges forming macroscopic objects. This makes it imperative to derive a theoretical formalism 
that is much simpler than the microscopic Maxwell—Lorentz electromagnetics and bypasses the 
unnecessarily detailed computation of the actual electromagnetic field. It turns out that doing this 
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is feasible by exploiting the two-layer structure of electromagnetics along with making hierarchal 
use of volume, time, and/or ensemble averaging. 
Indeed, in the words of Freeman Dyson [117], 


The modern view of the world that emerged from Maxwell’s theory is a world with two 
layers. The first layer, the layer of the fundamental constituents of the world, consists of 
fields satisfying simple linear equations. The second layer, the layer of the things that we 
can directly touch and measure, consists of mechanical stresses and energies and forces. 
The two layers are connected, because the quantities in the second layer are quadratic or bi- 
linear combinations of the quantities in the first layer. To calculate energies or stresses, you 
take the square of the electric field-strength or multiply one component of the field by an- 
other... The objects on the first layer, the objects that are truly fundamental, are abstrac- 
tions not directly accessible to our senses. The objects that we can feel and touch are on the 
second layer, and their behavior is only determined indirectly by the equations that operate 
on the first layer. 


Owing to this two-layer structure, the framework of the simplified theoretical formalism can be 
formulated as the following two-stage procedure: 


e first, define relevant optical observables as quadratic and bilinear forms in the electromagnet- 
ic field that (i) can be directly measured with suitably designed instruments, and/or (ii) quanti- 
fy the energy budget of a macroscopic object*; and 


e second, develop an efficient way to directly calculate appropriate averages of these observa- 
bles even if the detailed computation of the exact (microscopic) electromagnetic field itself is 
sacrificed. 


Indeed, the majority of applications do not require the knowledge of instantaneous (or quasi- 
instantaneous) local values of the optical observables but rather deal with averages taken over 
extended time intervals and/or finite (rather than infinitesimal) volume elements. Simple exam- 
ples of the experimental use of time averages are the exposure time of a camera, the integrating 
time of the rod cells in our eyes, and the signal integration over the sensitive face of a detector. 
Moreover, in many situations time averaging can be replaced by ensemble averaging, thereby 
resulting in further dramatic simplifications. 

In what follows, we will discuss how the two-layer structure of electromagnetics in combi- 
nation with appropriate averaging procedures yields an important effective-field approximation 
called macroscopic electromagnetics. This approximation is based on the introduction of a math- 
ematical entity called the macroscopic electromagnetic field and can be used to quantify time- 
harmonic electromagnetic scattering by a fixed macroscopic object. Although the macroscopic 
electromagnetic field in and of itself is not an actual physical field, it can yield suitably averaged 
optical observables directly, i.e., without the prior computation of the exact (i.e., microscopic) 
electromagnetic field. A straightforward generalization makes this approach applicable to quasi- 
monochromatic macroscopic fields and/or time-variable macroscopic objects. The resulting for- 
malism enables the computation of relevant time-averaged optical observables for a DRM by 
using analytical expressions and equations completely devoid of the actual electromagnetic field 
scattered by the medium. 


* Note that this definition of optical observables does not exclude quadratic and bilinear forms in 
the field that cannot be measured directly. The prime example of a bilinear form that often is not 
measurable is the Poynting vector [34]. 


All three types of averaging mentioned so far (i.e., volume, time, and ensemble averaging) 
have been used in various publications on electromagnetic scattering. Some of these publications 
may give the impression that different types of averaging may be used interchangeably or that a 
type of averaging can be selected almost at will. It is imperative to keep in mind however that 
each type of averaging has its own conditions of applicability and that indiscriminate use of any 
one of them can lead to physically meaningless results. Therefore, in what follows we will be 
very explicit in justifying the use of a specific type of averaging and explaining why an alterna- 
tive choice can be inappropriate. 


3. Macroscopic Maxwell equations 


War es ein Gott, der diese Zeichen schrieb (?) (Was it a God who wrote these signs (?)) 
[Ludwig Boltzmann [118], from Goethe’s Faust] 


From a long view of the history of mankind — seen from, say, ten thousand years from 
now — there can be little doubt that the most significant event of the 19th century will be 
judged as Maxwell’s discovery of the laws of electrodynamics. The American Civil War 
will pale into provincial insignificance in comparison with this important scientific 
event of the same decade. 

[Richard P. Feynman [119]] 


The macroscopic Maxwell equations (MMEs) were postulated by James Clerk Maxwell 
150 years ago [120] as the most fundamental laws of electromagnetics consistent with the totality 
of experimental data accumulated by that time. Maxwell’s ideas, summarized in his famous 
Treatise {|121], were picked up, systematized, and reworked mathematically by his immediate 
followers [122,123], most notably by Oliver Heaviside [124]. The subsequent notion that the 
MMEs must be a corollary of the more fundamental microscopic Maxwell—Lorentz equations 
[125] was put forth by Hendrik Lorentz in [112] and has been further developed by many authors 
[113—116,126—142]. 

In the framework of classical electromagnetics, the microscopic electromagnetic field is the 
only actual physical field which, in the vast majority of situations, is an extremely intricate func- 
tion of the position vector r and time ¢. The basic idea of macroscopic electromagnetics is that 
the detailed knowledge of the exceedingly complex dependence of the microscopic field on r and 
t is often not required in practice. Instead, this dependence is artificially simplified by averaging 
the microscopic field over either r or ¢, thereby yielding contrived macroscopic field vectors. It is 
imperative to recognize that these fictitious mathematical entities can only be useful to the extent 
to which they simplify the computation of macroscopic optical observables. In this respect mac- 
roscopic electromagnetics is the prime example of an effective-field approximation. 

The temporal variability of the microscopic electromagnetic field inside a macroscopic ob- 
ject is caused by the incessant microscopic movements of the constituent elementary charges, by 
macroscopic temporal changes of the object, and by time-harmonic oscillations and quasi- 
random fluctuations induced by the external sources. Throughout this Report, we will assume 
that field variations caused by macroscopic changes of the object occur much more slowly than 
those caused by the microscopic movements of the constituent elementary charges as well as 
much more slowly than the externally induced time-harmonic oscillations and, potentially, quasi- 
random fluctuations. Furthermore, we will assume that the quasi-random oscillations of the field 
occur much more slowly than its time-harmonic oscillations. A fundamental corollary of these 


assumptions is that monochromatic and quasi-monochromatic electromagnetic scattering by the 
slowly varying macroscopic object can be described at any moment in time by assuming that the 
object is fixed and solving the corresponding quasi-instantaneous boundary-value problem for 
the frequency-domain MMEs. 


3.1. Averaging over physically small volume elements 


The fundamental equations governing electromagnetic phenomena for point charges serv- 
ing as building blocks of a macroscopic material medium are the four microscopic Maxwell— 
Lorentz equations: 


V-e(r, 2) = ED (1) 
€0 
Vx e(r,f) + PUD — 9, (2) 
V - b(r,t) = 0, (3) 
de(r, t) ; 
Vx b(r, t) — €9Ho an a Mo jr, t), (4) 


where e and b are the microscopic electric and magnetic fields; 7 and j are the microscopic 
charge and current densities; €) and jg are the electric permittivity and the magnetic permeabil- 
ity of a vacuum; and 0 is a zero vector. Fundamentally, the microscopic fields are functions of r 
and ¢ only. This implies that the position vector and time are the only parameters over which 
e(r, 7) and b(r,t) can in principle be averaged. 

According to the first averaging approach dating back to Lorentz [112], the microscopic 
field is homogenized, at any moment in time, over “physically small” volume elements 6V cen- 
tered at r in order to smooth out drastic variations of e(r,t) and b(r,t) over interatomic dis- 
tances [114,126—128,133]: 


EB" (r,t) = (e(0, 0) 57, (5) 

BY (r,t) = (b(t, ))sv- (6) 

The well-known result of this approach is the system of the four MMEs: 

V-D" (r,t) =p" (9), (7) 
Vv 

VxE" (1,1) + BOD = 0, (8) 

V-B (r,t) =0, (9) 
OV 

Vx He (7, - PE = J’ (r,), (10) 


where H” (r,7¢) is the macroscopic magnetic intensity vector; p” (r, ¢) is the macroscopic elec- 


tric displacement vector; po” (r,t) and J wae t) are the macroscopic free charge density and 


current density, respectively. In the case of a macroscopically isotropic and time-dispersive ma- 
terial, the macroscopic field vectors entering the MMEs (7)—(10) are further related by the con- 


stitutive relations [116] 


t 
p*’ (r,1) = | dt’e™ (r,t — ’)E” (r, 1’), (11) 


: 1 
HH (r,t =| de’ —_—___ B*” (r, 7’), 12 
N= J dB) (12) 


t 
I” (r,t) = | dt’'o™” (r,t —1)E™ (r, 1). (13) 


One must recognize that this averaging procedure yields artificial field vectors formally 
satisfying the MMEs rather than an exact physical electromagnetic field. Furthermore, the prob- 
lem of actual practical significance is to compute macroscopic averages of specific optical ob- 
servables, including quantities describing electromagnetic energy budget. In other words, one 
needs volume averages of quadratic and bilinear forms in the electromagnetic field, such, for 
example, as the Poynting vector [143,144] 


s(r,t) = ee t) x b(r, ft). (14) 
0 


This implies that for a mathematical solution of the MMEs to be physically significant and prac- 
tically useful, it must enable the computation of relevant macroscopic optical observables, in- 
cluding the macroscopic Poynting vector 


Ss” (r, 2) = (s(t, D)ar = —- (e(r, 2) x BU, D) ar, (15) 


Ho 


by a simple substitution of the macroscopic field vectors for the microscopic ones, e.g., 


S$ (r,t) = -E(r,) xB (rd), (16) 
Lo 
Since the average of the vector product of two vectors in Eq. (15) is not necessarily equal to the 
vector product of the individual averages, Eq. (16) is highly nontrivial and by no means obvious. 
To the best of our knowledge, Eq. (16) has been derived rigorously only for structured pe- 
riodic nonmagnetic materials and only for the case of time-harmonic fields [141,142], while in 
all other situations (e.g., in the case of amorphous solids and liquids), it still has to be postulated. 
It must be recognized however that without Eq. (16) and similar formulas for other relevant sec- 
ond moments of the electromagnetic field the MMEs would largely lose their physical signifi- 
cance and become irrelevant or not helpful if one wishes to make useful predictions. 
Let us now assume that the macroscopic field vectors are monochromatic, while the medi- 
um is non-magnetic. It is convenient to represent the real field variables as real parts of the re- 
spective complex time-harmonic variables: 


E”’ (r,t) = Re[E* (r) exp(—ia)], (17) 
B” (r,t) = Re[B™ (r) exp(iaa)], (18) 
D” (r,t) = Re[D* (r) exp(—iat)], (19) 


HY (r,t) = Re[H™” (r) exp(—iar)], (20) 


p°" (r,t) = Re[p™” (r) exp(—ior)], (21) 
J’ (r,t) = Re[I™ (r) exp(-ior)], (22) 
where @ is the angular frequency, i = Ci > and “Re” stands for “the real part of”. It is then 


straightforward to show that the macroscopic field vectors satisfy the standard frequency-domain 
MMEs 


V-D" (r) =p (r), (23) 
Vx E (r)-iwB” (r) = 0, (24) 
V-B’ (r) =0, (25) 
Vx HA (r) + iad” (r) = J (r) (26) 
supplemented by the constitutive relations 

D” (r) = &" (r, @)E” (n), (27) 
H (r) = | ge (r), (28) 

Ho 
Jr) = 6" (r, aE” (1), (29) 


where the frequency-dependent electric permittivity € on (r,@) and conductivity gl (r,@) are, 
in general, complex valued: 


ze” (r,@) = if dt e% (r,t) exp(iot), (30) 
&” r,@) = J, dt o” (r,t) exp(iot). G31) 


According to Eq. (16), the near-instantaneous time-independent macroscopic Poynting vector is 
now given by the time average 


t+T/2 
s”(,.9) ==] 1 AEN (EF) AE) 
t-T 


= SRefE" (ry x [HA ON}, TT, (32) 
where the asterisk denotes a complex-conjugate value and 
i= a (33) 
o 


is the period of time-harmonic oscillations. 


3.2. Ensemble averaging 


The above methodology is intended to yield optical observables homogenized over physi- 
cally small volume elements at each moment in time. In the case of a high-frequency time- 
harmonic electromagnetic field, one could think of a different averaging approach intended to 
yield time-averaged optical observables at each point in space. Specifically, the microscopic 
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electric and magnetic fields as well as the microscopic charge and current densities are factorized 
according to 


e(r,t) = Re[e(r, t) exp(—iaf)], (34) 
b(r, 2) = Re[b(r, t) exp(-iar)], (35) 
n(r,t) = Re[7(r, t) exp(—ia@t)], (36) 
j(r, 0) = Ref j(r, 0) exp(—iar)| (37) 


[131,133]. The dependence of the complex amplitudes e(r, f), b(r, 0), n(r,t), and ja.o on 
time is implicit in that it is caused by relatively slow random movements of the microscopic 
charges occurring independently of the rapid oscillatory motions described by the time-harmonic 
factor exp(—ia@t). This means that at any moment in time, the complex amplitudes satisfy the 


frequency-domain microscopic Maxwell—Lorentz equations 


Vite pe 2, (38) 
€0 
V x @(r, t) — iwb(r, t) = 0, (39) 
V-b(r, ft) = 0, (40) 
V x b(r,f) + i@eg up € (r,t) = Lo (r,1) (41) 
provided that 
de(r,t)| sewed! (42) 
ot 
Pred < ofb(r,o| (43) 


Note that whether the inequalities (42) and (43) are satisfied can be expected to depend on sever- 
al factors, including the angular frequency, the material in question, and the material tempera- 
ture. 

Averaging the fields (34) and (35) over time is meaningless since the rapidly oscillating 
factor exp(—i@t) causes both averages to vanish: 


lp etT/ , os 
=| dt exp(-iot) = 0. (44) 
T J t-T/2 T > T. 


° 


However, the vector product of the complex electric field and the complex conjugate of the mag- 
netic field varies with time much more slowly since the factors exp(—iw@t) and [exp(—iat)]" 
cancel each other. Therefore, averaging e(r,t) x b(r,t) over a period of time 7 much longer than 
T, but much shorter than the typical temporal scale T’ of variability of the complex amplitudes 
€(r,t) and b(r,t) is meaningful and yields a quasi-instantaneous Poynting vector slowly vary- 
ing in time: 

i. 1 


t+T/2 : P P 
((s(r,t))) = ae dt’e(r, t’) x b(r, 1’) 


te —— Refe(r,f)x[b(r, AT}, ).K TT’, (45) 
Ho 


where the symbol ((---)) hereinafter denotes averaging over a sufficiently long time interval, the 
actual length of the time interval being clear from the context. The time-independent macroscop- 
ic Poynting vector is now defined as the average over a time interval T much longer than 7” : 


1 1 ie ~ , 1 n7* , 
((S(r, £))) = ——— Re dfe(r,f)x[bar, OH], TT’. (46) 
2Uo T t-T/2 
The computation of ((S(r, f))) is usually simplified by assuming ergodicity of the ensemble 
of elementary charges (see, e.g., Section 10.4 of Ref. [34]°) and replacing the temporal average 
in Eq. (46) by the statistical average over the ensemble y of configurations of all the microscop- 
ic charges: 


Sr.) = Re J dy Bey) x1BU0, WAT PLY) (47) 
0 


where p(y) is a suitable time-independent probability density function. Similar expressions can 
be written for macroscopic versions of other quadratic and bilinear forms in the microscopic 
electromagnetic field. 

In general, the computation of the ensemble average (47) is still quite involved. A major 
simplification is based on the as yet unproven assumption according to which 


(S(r, 8) = ——RefE" (r) x [BY (NT) (48) 
29 


(and similarly for other second moments), where all materials are assumed to be non-magnetic 
and the corresponding macroscopic complex field vectors are given by the individual time aver- 
ages replaced by ensemble averages: 


- 1 t+T/2 _ a 

EY (r) = — | dt’ €(r,1’) = | dye(r,y)py), TT", (49) 
T J t-T/2 

S Lp ttT/2 7 a 

BY(n=—[ dr bar) =[dybewpw, T>T (50) 
T 4 t-T/2 


It can then be shown [131] that the macroscopic field vectors satisfy the standard frequency- 
domain MMEs 


V-D’(r) = p(n), (51) 
V x EY (r)-i@B” (r) = 0, (52) 
V- BY (r) = 0, (53) 
V x HY’ (r) + i@D” (r) = JY (r) (54) 


> Instructive discussions of the ergodic hypothesis as a basic underlying principle of classical and 
quantum statistical mechanics and kinetic theory can be found in [145-149]. 
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supplemented by the constitutive relations 


DY (r) = &"(r, o)E" (r), (55) 

H” (r) = | BY, (56) 
Ho 

J” (r) = 6 (r, @)E" (r), (57) 


where, again, the frequency-dependent macroscopic electric permittivity €” (r,@) and conduc- 
tivity c” (r,@) are, in general, complex valued. 

We have seen in the preceding subsection that equations having the same mathematical 
structure as Eqs. (51)-(57) can also be obtained using the volume-averaging approach. As a con- 
sequence, it is often believed that a formal re-multiplication of the vectors EY (r) and BY (r) 
given by Eqs. (49) and (50) by the time-harmonic factor exp(—i@t) yields the ensemble- 
averaged time-dependent electromagnetic field. This belief is questionable since ensemble aver- 
aging is not a primordial physical concept and can only be used as a substitute for time averaging 
(see, e.g., [146] and pages 1-6 of [150]). Averaging the right-hand sides of Eqs. (34) and (35) 
over time yields a zero result provided that T, < 7’. Therefore, the ensemble-averaged time- 
dependent electromagnetic field must also be zero. The reader should recall that the vectors 
E’ (r) and BY (r) are obtained by: 


e artificially omitting the time-harmonic factor exp(—ia@t) in Eqs. (34) and (35); 


e taking the time average of the remaining factors; and 


e replacing this time average by the ensemble average based on the assumption of ergodicity 
according to Eqs. (49) and (50). 


It is thus obvious that the subsequent re-multiplication of E’(r) and BY(r) by exp(—iat) 
yields quantities of questionable veracity rather than actual time or ensemble averages. 


3.3. Further discussion 


The ensemble averaging approach described in Subsection 3.2 bypasses the introduction of 
time-dependent macroscopic vector fields altogether, whereas the volume averaging approach 
outlined in Subsection 3.1 does yield macroscopic field vectors explicitly dependent on time as 
well as on coordinates. As a consequence, certain solutions of the time-domain MMEs (7)-(10) 
do describe vector waves unfolding in space and time, which may seem to provide substance to 
the widespread belief that E* (r,t) and H® (r,t) represent an actual physical field in a “ho- 
mogenized macroscopic medium” rather than a purely mathematical entity. It is imperative to 
recognize however that (i) any type of averaging is a purely human intervention resulting in an 
artificial rather than an actual physical field; (ii) the computation of the macroscopic field vectors 
E” (r,t) and H® (r,t) (or EY (r) and BY (r)) is almost never an end in itself; and (iii) irre- 
spective of the averaging approach chosen, the ultimate purpose of macroscopic electromagnet- 
ics is the computation of time- and/or volume-averaged optical observables and the average elec- 
tromagnetic energy budget. The MMEs are useful and meaningful only to the extent to which 
they help achieve this objective by eliminating the need to solve explicitly the microscopic 
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Maxwell—Lorentz equations. Despite substantial recent progress in the microphysical justifica- 
tion of the MMEs, this problem still awaits a definitive solution. 

Although the formal mathematical structure of Eqs. (23)-(29) and (51)-(57) is the same, 
their solutions can, in principle, be different. Indeed, the specific procedures for the computation 
of the corresponding macroscopic electric permittivities and conductivities are not necessarily 
the same, and it is far from being obvious that €°” (r,@) = €” (r,@) and 6°” (r,@) = 6” (r, 0). 
In what follows, we will usually omit the superscripts 6V and wy for the sake of brevity, but it 
should be kept in mind that the actual values of €(r,@) and o(r,@) may depend on the averag- 
ing approach chosen. 

The relative merits of either homogenizing optical observables over physically small vol- 
ume elements or averaging them over a time interval J >>T’ remain somewhat unclear. The two 
averaging strategies are likely to yield similar results if they are applied to the calculation of the 
time-averaged radiation budget of a macroscopic volume with dimensions greatly exceeding the 
wavelength, provided that @°"(r,w) =@”(r,@) and &°"(r,w) = G”(r,@). However, the 
modeling of the interaction of the electromagnetic field with a photodetector may depend on 
whether one uses the macroscopic field vectors homogenized over physically small volume ele- 
ments at a given moment in time or those averaged over time at a given point in space. This issue 
obviously needs to be further clarified. 

On the more fundamental level, the MMEs must be derived in the framework of the QED 
by quantizing the microscopic electromagnetic field. There has been notable progress in this di- 
rection [151—154], but definitive studies are still needed. 


4. Monochromatic and quasi-monochromatic scattering by a fixed macroscopic object 


Consistent with the preceding discussion, the foundation of the theory of electromagnetic 
scattering by a DRM can be built of the following four major building blocks: 


the theory of monochromatic scattering by a fixed finite object; 


the theory of quasi-monochromatic scattering by a fixed finite object; 


the theory of monochromatic scattering by a randomly changing object; and 


e the theory of quasi-monochromatic scattering by a randomly changing object. 


The main purpose of this section is to give an explicit formulation of the electromagnetic scatter- 
ing problem for a fixed object in maximally general terms and discuss its immediate implica- 
tions. We start with monochromatic scattering and then, in Subsection 4.12, generalize the for- 
malism to encompass the case of quasi-monochromatic radiation. Monochromatic and quasi- 
monochromatic scattering by a time-variable object such as a DRM will be considered in the fol- 
lowing section. 

Let us define the characteristic length / according to 


dK«Il<), (58) 
where d is the typical distance between a molecule and its closest neighbors and / is the wave- 


length. A major premise of the previous discussion is that frequency-domain macroscopic elec- 
tromagnetics can be used to 


e compute the value of any second moment of the microscopic electromagnetic field homoge- 
nized over physically small volume elements with dimensions of the order of / and averaged 
over a time interval 7 >>7, (hereinafter averaging strategy 1; Subsection 3.1), or 


e compute the value of any second moment of the microscopic electromagnetic field averaged 
over a time interval T>>7” at any fixed point in space (hereinafter averaging strategy 2; Sub- 
section 3.2). 


This is done by first solving the frequency-domain MMEs 


V- D(r) = P(r), (59) 
V x E(r) —i@B(r) = 0, (60) 
V - Br) = 0, (61) 
V x H(r) + i@D(r) = J(r) (62) 
supplemented by the constitutive relations 

D(r) = @(r, @)E(r), (63) 
f(r) = Bin, (64) 

Ho 
J(r) = G(r, @)E(r) (65) 


and then substituting the resulting macroscopic field vectors E(r) and B(r) in the appropriately 
modified formula for the requisite second moment of the microscopic field. Typically the modi- 
fication amounts to applying the operation “(1/2) Re”. For example, the microscopic Poynting 


vector is defined by Eq. (14), while its time-independent average macroscopic counterpart is giv- 
en by 


Sir) = 5 Ret Btn) x [A]"}, (66) 


where S(r) stands for either (S°” (r,4)) or ((S(r, £))). According to [155], the frequency-domain 
macroscopic formalism can be expected to work well as long as the smallest homogeneous ele- 
ment of the scattering object exceeds ~50A. And even for smaller elements, it may produce 
meaningful results if combined with empirical corrections for €(r,@) and o(r,@) [156]. 
The three basic ingredients of some phenomenological approaches to electromagnetic scat- 
tering by a DRM have been: 
(i) the visual perception of the DRM as being assembled of separate “building blocks” in the 
form of particles; 
(ii) the presumed knowledge of how to compute specific optical observables for each individu- 
al building block in the absence of all the other blocks; and 
(iii) the belief that the optical observables for the assembly of the building blocks can somehow 
be expressed in terms of the optical observables computed for the separate building blocks. 


The latter belief has usually been justified by verbal “simple physical considerations” and ac- 
cepted as needing no rigorous mathematical proof. However, this approach is generally incorrect 
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since unlike the human eye, the electromagnetic field perceives the DRM in its entirety rather 
than one particle at a time. Therefore, any first-principles approach to electromagnetic scattering 
by a DRM must originate in the explicit formulation of the MMEs and appropriate boundary 
conditions for the entire DRM rather than in the set of separate formulations for the individual 
particles. 


4.1. The standard scattering problem 


Consider a fixed finite scattering object imbedded in an infinite medium; the latter is as- 
sumed to be homogeneous, linear, isotropic, nonmagnetic, and nonabsorbing. In general, the ob- 
ject is a cluster consisting of a finite number N >1 of separated or touching distinct compo- 
nents. It occupies a finite interior region V;y7 given by 


N 
Vny = U FV, (67) 
i=1 


where V, is the interior volume of the ith component (Fig. 6). The object is surrounded by the 


infinite exterior region Vgy_ defined such that Viz U Vexp = K°, where KR? denotes the entire 
three-dimensional space. It is further assumed that the interior volume Vjy7 is filled with iso- 


tropic, linear, nonmagnetic, and possibly inhomogeneous material. Point O serves as both the 
common origin of all position vectors and the origin of the laboratory coordinate system. 

Unlike the general microscopic Maxwell—Lorentz equations (1)—(4), the four frequency- 
domain MMEs (59)-(62) are not independent since Eqs. (59) and (61) follow from Eqs. (60) and 
(62) [34]. This allows one to consider only the Maxwell curl equations, re-written as 


V x E(r) =i@pyH(r) 


Z U i. 68 
| ae 8) 


V x E(r) = iopyH(r) 


V x H(r) = ee re Vr, (69) 


where ¢; is the real-valued electric permittivity of the infinite host medium and 


S5(t,@) = B_(r,0) + j 220) (70) 


is the (potentially co-ordinate-dependent) so-called complex permittivity of the scattering object. 
The corresponding boundary conditions read: 
fi x [E,(r) — E,(r)] = 0 


1 
i x [fi (r) — A (r)] = i re Sint? om 


where the subscripts 1 and 2 correspond to the exterior and interior sides of the boundary Sy 
of the object, respectively, n is the local outward normal to S;y;, and Sj,y7 is the union of the 
closed surfaces of the N components of the object: 


N 
Sint = U S;. (72) 
i=l 


Let us now assume that the complex amplitudes E(r) and H(r) everywhere in %* can be 
written as a superposition of the complex amplitudes of a plane-wave “incident field” (super- 


script “inc”) propagating in the direction of the unit vector i'™° and those of a “scattered field” 
(superscript “sca’”): 


E(r) = B™(r) +E (, =) 
H(r) = H'°(r) + HS (r), = 
where 

E"<(r) = Ex exp(ikn'™ . r), (75) 
H'(r) _ Hi’ exp(ikn™ . r) = ae nize x Eire exp(ikn'™® . r), (76) 

V Lo 
and 
ky =@Vé\Lo me 


is the wave number in the exterior region Vgx;. Furthermore, we require the scattered field am- 
plitudes to satisfy the following asymptotic condition at infinity: 


lim {4/9 r x HS“ (r) + role, ES? (r)} = 0, (78) 
r— 


where r =|r | is the distance from O to the observation point (Fig. 6). The limit (78) is tradi- 
tionally called the Silver—Miiller radiation condition at infinity [157,158] and holds uniformly 
over all outgoing directions r = r/r. 

The combination of the Maxwell curl equations (68) and (69), the boundary conditions 
(71), the decomposition (73)-(76), and the asymptotic condition (78) represents the so-called 
standard electromagnetic scattering problem for plane-wave illumination. 

The mathematical decomposition (73)-(76) of the “total” macroscopic frequency-domain 
field vectors makes it clear that the scattered field is defined as the difference between the total 
fields corresponding to the situations when the object is and is not present. This is consistent 
with the point of view that the incident field is not transformed into or replaced by the scattered 
field. In other words, the physical cause of frequency-domain scattering by the object is not the 
incident field, but rather the very presence of an object with optical properties different from that 
of the exterior medium [159,160]. 

Note that to the best of our knowledge, the boundary conditions (71) have not been derived 
explicitly from the microscopic Maxwell—Lorentz equations. Ad hoc ways of introducing these 
conditions in the framework of macroscopic electromagnetics are discussed in Section 2.8 of 
[161]. 

Since the first relations in Eqs. (68) and (69) yield the magnetic field vector provided that 
the electric field vector is known everywhere, the solution of the standard scattering problem is 
often sought in terms of only the electric field vector. 


4.2. Existence and uniqueness of solution of the standard problem 


The statement of the standard scattering problem would be of little practical use if this 
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problem had no solution and/or if the solution was not unique. Fortunately, both the existence 
and the uniqueness of solution have recently been demonstrated for particles with smooth surfac- 
es (see [162,163] as well as Section 9.1 of [115]). Certain results for particles with edges do ex- 
ist, but this case is fundamentally more difficult since the formulation of the boundary condition 
becomes ambiguous unless appropriately modified (see Chapter 9.2 of [115] and the discussion 
in Subsection 4.3). 


4.3. Volume integral equation 


Although the standard scattering problem is formulated for the incident field in the form of a 
plane electromagnetic wave, the range of its actual applicability is much wider. Indeed, the line- 
arity of both the MMEs, the boundary conditions, and the radiation condition at infinity implies 
that solving the standard problem yields the solution of a more general scattering problem as 
long as the corresponding incident electromagnetic field can be expanded in plane electromag- 
netic waves. This profound fact becomes especially evident if we consider a mathematically 
equivalent formulation [162,163] of electromagnetic scattering in terms of the so-called volume 
integral equation (VIE) [34,164] (see also [165]): 


E(r) = E™(r) + 7 | @r’G(r, r’) - Er’) [m7 (r’) — 1] 
Vint 
exp(ik,|r — r’) 


_ E'™ (r) de ai +4 ra ® v) s | d3r’ E(r’) m’(r’) _ 1], re KR, (79) 
1 


Vint 4a Ir = 


where 


beat [7 sive v) expGiki|r — rl) (80) 
ke An |r —r"| 


is the free-space dyadic Green’s function, 


Mss? (r’, @) 


(81) 
&| 

is the (complex) refractive index of the object’s interior relative to that of the host exterior medi- 
um, / is the identity dyadic, and @ is the dyadic product sign. 

The second equality of Eq. (79) is a mathematically rigorous expression which has been 
used in [162,163,166] to deduce several useful corollaries. By contrast, the first equality is a 
shorter expression, but contains a strong singularity (strictly speaking, a non-integrable one) 
when re V;yr. Then the integration must be carried in the following specific principal-value 
sense to ensure that it is equivalent to the rigorous expression [167,168]: 

3A , , : 3A , , i47 
| drG(r,r)F(r) = lim d'rG(r,r)F(r) — — F(r), (82) 
Vint Yo 30 J Vint \Y 3 

where V» is a spherical exclusion volume around r. In what follows, we use the compact version 
of Eq. (79) and similar ones, but always assume that it implies the abbreviation (82). 

The VIE incorporates the boundary and radiation conditions and expresses the total field 


everywhere in space in terms of the total internal field. It can even be considered to be more in- 
clusive since it is well behaved (has a unique solution) even for particles with sharp edges [166]. 
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Therefore, in the following we do not impose any limitations on the object’s boundary and as- 
sume, in a somewhat ad hoc fashion, that in the presence of sharp edges the scattering problem is 
formulated through its VIE representation. 

Owing to its specific mathematical form, the VIE serves as the very embodiment of the 
concept of frequency-domain electromagnetic scattering [159,160]. Indeed, it shows that if the 
scattering object is absent (m(r’) = 1), then the total field is identically equal to the incident 


field. The presence of the object (m(r’) # 1) changes the total field, thereby allowing the defini- 


tion of the scattered field as the difference between the total fields in the presence and in the ab- 
sence of the object. Furthermore, the VIE implies that the incident field is not modified by the 
presence of the object and, thus, is not transformed into the scattered field. 

The linearity of the VIE suggests that it should be convenient in many cases to express, pure- 
ly mathematically, the scattered electric field in terms of the incident field: 


ES (r) = | 


where 7 is the so-called dyadic transition operator of the scattering object. Eqs. (79) and (83) 


imply the following integral equation for 7 traditionally called the Lippmann—Schwinger equa- 
tion (cf. [10,169,170]): 


@r’G(r, 1) | Pr’'T (rr) E*(r"), re 8, (83) 


Vint Vint 


T(r,r’) = kp [m7 (r) — 1] 6(r - rv’) 


+ k?[m?(r) - 1] | @r’G(r,r’)-T(r’,r’), or,r € Vint, (84) 
Vint 


where 5(r) is the three-dimensional delta function. Note that 7 is independent of the electro- 
magnetic field and is defined only by the spatial distribution of the relative refractive index 
throughout V;,7. As such, it can be viewed as a complete “optical identifier” of the scattering 


object. 


4.4. Scattering in the far zone of the entire object 


The spatial distribution of E(r) and H(r) inside the scattering object as well as in its im- 
mediate vicinity can be quite complex. However, there is a drastic simplification as the distance 
from the object increases since, irrespective of the specific nature of the object, the scattered 
field ultimately becomes a spherical outgoing electromagnetic wave. Indeed, a key property of 
the dyadic Green’s function is the asymptotic behavior 


exp(ik,r) 


exp(—ifir - r’), (85) 
Anr 


G(r,r’) > (1 -f @f) 
roo 


where, as before, r = r/r. As a consequence, placing the origin of the laboratory coordinate sys- 


tem O at the geometrical center of the scattering object, as shown in Fig. 7, and substituting Eqs. 
(75) and (85) in Eq. (83) yields [34] 


Es? (r) = exp(ikr) E¥? (n°) ~ exp(ikir) A(n, ni’) . EC ns“. Ej? (n°) =) (86) 
r r 


roo 


Asca inc 


Here, n°“ = F is a unit vector in the scattering direction and Ans, ni") is the scattering dyad- 
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ic having the dimension of length and possessing the properties 
n@ AG”, ni”) = —_ Aas, n™) 2 nie = 0. (87) 


The explicit expression for the scattering dyadic in terms of the dyadic transition operator is as 
follows: 


An, n'”’) _ =f. n°? ® n°“). I, d’r’exp(—ikn*" . r’) 
An 
x| Br’ T(r’, r’)exp(ikin'™ - r”)- (7 —n'™ @ ni). (88) 
Vint 


Eq. (86) implies that the electric and magnetic field vectors of the scattered electromagnetic field 
vibrate in the plane perpendicular to the propagation direction and decay inversely with distance 
from the object. 

The formal mathematical conditions of applicability of Eq. (86) are as follows: 


K(r-a)>1, (89) 

r>a, (90) 
, 

r>> nee, (91) 


where a is the radius of the smallest circumscribing sphere of the entire scattering object cen- 
tered at O. The physical meaning of these inequalities is discussed in [34,171]. 

The main attraction of the far-zone approximation is that the entire object is implicitly 
treated as a point source of scattered radiation, while the scattered field is reduced to a simple 
outgoing spherical wave. Furthermore, Eq. (87) implies that out of the nine components of the 
scattering dyadic in spherical coordinates centered at the origin, only four are independent. It is 
thus convenient to introduce a 2x2 amplitude scattering matrix S having the dimension of 
length and expressing the 9- and @-components of the scattered spherical wave in the 9- and 


gy- components of the incident plane wave: 


Ee" (rns) exp(ikir) a. : Eat 
S(n a nize 92 
oe r (n ) Eye | ( ) 


Here, 6 € [0, z] is the polar (zenith) angle measured from the positive z-axis and @ ¢€ [0, 27) is 
the azimuth angle measured from the positive x-axis in the clockwise direction when looking in 
the direction of the positive z-axis. 

A fundamental property of the scattering dyadic and the amplitude scattering matrix is the 
following symmetry with respect to reversing and interchanging the incidence and scattering di- 
rections [172]: 


A(-n"’, _ n°) = [An nize ee (93a) 
S(-i™, _ n’) = = : os ve 7 Sa ia’ n inc | (93b) 
Si2(f ) Soo (i ) 


where T stands for “transposed.” Eqs. (93a) and (93b) are traditionally called reciprocity rela- 
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tions. 


4.5. Well-collimated radiometers 


By solving the MMEs either analytically or numerically, one can model a wide range of 
optical observables, including those that can be measured with actual optical instruments. Some 
of these instruments are expressly intended for near-field applications [173], while some can 
measure both near- and far-field observables. As explained in [100,174], the overwhelming ma- 
jority of laboratory, in situ, and remote-sensing instruments measuring specific manifestations of 
electromagnetic energy transport in particulate media belong to the category of well-collimated 
radiometers (WCRs). Depending on the measurement setting, these instruments can work in the 
near as well as in the far zone of the particulate scattering object, but in either case it is assumed 
that the total electromagnetic field at the observation point is a superposition of plane or near- 
plane wavefronts. 

The principal functional elements of a WCR are the objective and relay lenses, the pinhole 
diaphragm, and the photoelectric detector, as shown schematically in Fig. 8a. The physical na- 
ture of the measurement afforded by the WCR can be illustrated by considering the response of 
the instrument to the field formed by superposing two plane electromagnetic waves propagating 
in the directions of the unit vectors q, and q,, respectively. The effect of the objective lens on 
the total field is a superposition of its effects on each plane-wave component. According to the 
paraxial approximation (see Section 5.1 of [175]), in the near zone of the objective lens either 
plane wavefront is transformed into a converging spherical wavefront (Fig. 8b) with its focal 
point located in the plane of the diaphragm. The first spherical wavefront passes freely through 
the pinhole, is converted back into a plane wavefront by the relay lens, and impinges on the sen- 
sitive surface of the photodetector, thereby defining the signal generated by the WCR. The sec- 
ond spherical wavefront becomes extinguished by the diaphragm and never reaches the photode- 
tector. Thus the combination {objective lens, diaphragm} serves to filter out only plane (or near- 
plane) wavefronts propagating in directions very close to the optical axis of the WCR and falling 
within its small acceptance solid angle 

a 
AQ =, (94) 
a 
where d is the diameter of the pinhole and fis the focal length of the objective lens. 

Typically a photodetector reacts only to the intensity of the beam impinging on its sensitive 
surface. However, by inserting special optical elements (such as polarizers and retarders) be- 
tween the relay lens and the detector in Fig. 8a, it is possible to modify the resulting beam im- 
pinging on the detector in such a way that its new intensity contains information about the polar- 
ization state of the original superposition of the plane or near-plane wavefronts filtered out by the 
{objective lens, diaphragm} combination. The result is a photopolarimetric WCR. 

Despite having quite different appearances, the one natural and six manmade devices in 
Fig. 9 have the same basic physical function: they filter out electromagnetic wavefronts rather 
than electromagnetic energy currents. In a radio telescope (Fig. 9b) or a reflecting optical tele- 
scope (Fig. 9c), the functional role of the objective lens is played by the radio antenna or the 
primary mirror, respectively. In the final analysis, however, all these devices are WCRs, possibly 
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with an added panoramic (or imaging) capability’. The basic functionality of a WCR makes it 
quite useful in cases when the total electromagnetic field at the observation point can naturally 
be represented as a superposition of plane or locally near-plane wavefronts, the electromagnetic 
field in the far zone of a finite object being a prime example. Since the measurement enabled by 
a WCR is well defined in terms of basic concepts of light—matter interactions, it should be ame- 
nable to theoretical modeling. This explains why the combination of a WCR and an appropriate 
theoretical analysis tool often serves as an efficient means of optical characterization. 


4.6. Far-zone optical observables 


The formalism embodied by Eqs. (75) and (92) helps illustrate how to define specific far- 
field optical observables measurable with WCRs. The main results of the following analysis will 
be straightforward consequences of the total field in the far zone being a superposition of two 
transverse wavefronts, i.e., the incident plane-wave field and the scattered spherical-wave field. 

As we have already mentioned, interposing one or more optical elements between the relay 
lens and the photodetector of a WCR can enable it to measure the power corresponding to partic- 
ular polarization components of the superposition of plane or near-plane wavefront passed by the 
{objective lens, diaphragm} angular filter. Similarly, interposing one or more such optical ele- 
ments before the scattering object, we can generate the incident field with a specific state of po- 
larization. Repeating the measurement for a number of different combinations and/or orientations 
of these optical elements enables us to determine the specific mathematical relationship between 
a complete set of polarization characteristics of the incident field and that of the field impinging 
on the objective lens of a WCR. This relationship is traditionally formulated in terms of 4- 
element columns formed by the Stokes parameters and 4x4 so-called phase and extinction ma- 
trices. 

According to the preceding discussion, a WCR that is not facing the incident wave and is 
not centered at the scattering object will not generate any signal. Therefore, let us first consider 
the situation when the instrument has its optical axis centered at the object in the scattering direc- 


tion away from the incidence direction, i.e., F#m" (WCR 1 in Fig. 10). It is clear that in this 
case the instrument filters out only the quasi-plane part of the outgoing spherical wave cut out by 
its objective lens, as shown schematically by the dashed curve. Therefore, the average polariza- 
tion response of WCR | per unit time can be expressed in terms of the so-called Stokes column 
vector of the scattered wave as follows: 


Signal 1 = S,1°“(rn**), (95) 


where the overbar has the same meaning as in Eq. (66), r is the distance from the scattering ob- 
ject to WCR 1, n“" =r, and S,, is the surface area of the objective lens. Recalling the defini- 


tion of the real-valued Stokes parameters of a transverse electromagnetic wave [20,34] and Eq. 
(92), we have 


® In this case each pixel of a charge-coupled device or each photoreceptor cell of the retina has a 
dual role of the diaphragm and the detector of electromagnetic energy flow [34]. 
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(96) 


Analogously, the polarization state of the plane incident wave (75) can be characterized in terms 
of the Stokes column vector 


pre ae (EY a Ee (Ei Mo) 

. inc Eee Ejmey* = Ei Einey* 
uy" 2 Ho _ Ear (Ev) = i (Ey 
yr [Ege Gy 22 Bee Euey: | 


The es transformation law reads: 
18 (r ac) = nH, cor (98) 
r 


where Z(n°, a'™°) is the 4x 4 Stokes phase matrix with real-valued elements given by the fol- 
lowing wine and bilinear combinations of the elements of the amplitude scattering matrix 


S(a’*, A?) [20,25,34]: 


21 = 50 Su P +1 Si P +1 So P +1 Soe P), (99) 
Z\2 = sll Sir P =| Siz P +] Sar P =| S22 7), (100) 
Z)3 = —Re(S\Sj2 + SxS), (101) 
Zy4 = —Im(S\)S}9 — SxS91), (102) 
Zn = 5S P +1 S2 P -| Sar P -1 S22 P), (103) 
Zm = 5 Su P -1 SP - 18a P +1 S22 P), (104) 
Z3 = Re(S})S}2 — S39531), (105) 
Z4 = —Im(S})S}2 + Sy7S51), (106) 
Z31 = —Re(S1183 + So95}2), (107) 
Z39 = — Re(S\)S31 — Sy2Sj2), (108) 
Z33 = Re(S\S39 + S12S31), (109) 
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Z34 = Im(S,,S 7) + S151), (110) 


Z41 = —Im(Sy)S}, + SoS), (111) 
Zay = — Im(S31S}1 — S225}9), (112) 
Z43 = Im(Sy7S}1 — $1253), (113) 
Z44 = Re(Sp9S{1 — Sy2S91). (114) 


The Stokes phase matrix has the dimension of area. 
Let us now consider a polarimetric WCR with its axis centered at the scattering object in 


the exact forward-scattering direction f = n'™°, i.e., WCR 2 in Fig. 10. Now the {objective lens, 


diaphragm} angular filter of the instrument accepts both the incident plane wave and the part of 
the outgoing spherical wave propagating in the forward-scattering direction and cut out by the 
objective lens. As a consequence, we can define the Stokes column vector of the total field for 


propagation directions f very close to a'™° 
Eg(rP Eg (rF)] + Eg (PEG?) 
Eg(r? Eg?) — Eg (FE QV? 
I(rt’) = Ta alr zy alr Zl ol” J ol Q wb be AQ), (115) 
2 Vio] — Eo PEP — EP LEG"? 
Li{Eg (PE g(r FY — Eg (PE (rP YT 3 | 
where AQ, is the acceptance solid angle of WCR 2 and the total electric field is given by 


E(r'f’) = E'°(r'f’) + ES (r'P). (116) 


war 


Integrating the elements of I(7'F’) over the objective lens of WCR 2 yields the following expres- 
sion for the average recorded polarized signal per unit time [20,25,34]: 


Signal 2 _ Sirs Ki ye + Sot Zn, Or”. (117) 
r 


ae a nan 


where Z(n ) is the forward-scattering Stokes phase matrix and K(f'"°) is the real 4x 4 


Stokes ence matrix. The elements of the latter are linear combinations of the elements of 


the forward-scattering amplitude matrix S(a'"°, nl): 


2a : 
Ky = eo + S59), J= 1,..., 4, (118) 
20 
Kyy = Ky, = 7 en — S59); (119) 
1 
20 
Ky3 = K3) = ae nO + So1), (120) 
1 
20 
Ky4q = Ka, = PCr — Siz), (121) 


1 


26 


20 
Ky3 = —K32 = 4, a — S12); (122) 
1 
20 
Koq = —Kay = a + S51), (123) 
1 


20 
K34 = -K43 = a, See _ Si1)> (124) 
1 


where “Im” stands for “imaginary part of’. Like the phase matrix, the extinction matrix has the 
dimension of area. 

Eq. (117) is the most general form of the so-called optical theorem. It demonstrates that the 
presence of the scattering object not only changes the total power of the electromagnetic radia- 
tion recorded by the WCR facing the incident wave (WCR 2 in Fig. 10), but also can change its 
state of polarization. The latter phenomenon is called dichroism and is caused by different atten- 
uation rates for different polarization components of the incident wave in the case of an object 
lacking perfect spherical symmetry. Moving WCR 2 sufficiently far from the scattering object 
can render the contribution of the third term on the right-hand side of Eq. (117) negligibly small, 


Signal2 = S,1'™ -K(a'™)yr", (125) 
r3© 
and thereby make the extinction matrix a directly measurable quantity. 


Among the most general properties of the phase and extinction matrices [18,20,178] are the 
inequalities 


CS fi, (Zyl Ai fa hhc): (126) 
the reciprocity relations 

Z(—-n"™’,-n) = A3[Z(a*, H™)]" Ag (127) 
and 


K(-a™) = As[KGi"™)]" As; 
the backscattering theorem 
Z\\(-M, N) — Z2(—M, M) + Z33(—M, Nn) — Z44(—n, ) = 0; (128) 
and the symmetry relation 
K, ; (a’"*) Ki (a) _ Kis (a’"*) Kia (a"*) 
Ky (a) Kn (A) Ko3 (A) — Koa (a) 
_ K; 1 (n'"*) Ky (n°) K33 (a) Ky (a) 
Kay (n'"* ) _ Ka (n'"* ) Ka; (n'"* ) Kay (n'"* ) 


K(-f'™) = (129) 


where 
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1 0 0 
0 1 0 

A, = Al =A;' = (130) 
0 =-1 0 
00 0 1 


The properties (127)—(130) follow directly from the reciprocity relation (93) combined with Eqs. 
(99)-(114) and (118)-(124). 


4.7. Energy budget 


In the preceding subsection we explained how to quantify the electromagnetic response of a 
far-field WCR. In this subsection we discuss the theoretical solution of the energy-budget prob- 
lem for an arbitrary volume V enclosing the entire scattering object (Fig. 11a). Since the host 
medium is assumed to be nonabsorbing, the net average rate at which the electromagnetic energy 
crosses the closed boundary S of the volume is always nonnegative and is equal to the power ab- 
sorbed by the object: 


vs - -| @r S(r) «A, (131) 
Ss 

where, as before, 

S(r) = + Re{Bir) x AO} (132) 


is the average macroscopic Poynting vector and n is the unit vector in the direction of the local 


outward normal to S. According to Eqs. (73) and (74), ws can be written as a combination of 
three terms: 


wes — prime _ sea, pest (133) 
where 

We = ~zRe [ rE" @) xf OT} a (134) 
ys _ = Re I. ar (ES? (p) x [H**(r)]"} ‘ n, (135) 
Ww = : Re I. @r (pine (r) x (Hs (r)J" + Es? (r) x (Hi"’ (r)]"} . n. (136) 


It can easily be seen that W'° vanishes identically because {E"°(r) x [H'™(r)]"}/2 is a constant 


vector independent of r, which is a trivial consequence of the surrounding medium being loss- 
less. Therefore, the absorption rate is equal to the difference between the energy extinction rate 
and the energy scattering rate: 


ws = we _ ws. (137) 


Again, one can exploit the assumption that the infinite host medium surrounding the object is 


nonabsorbing to show that the values of W°' and W ‘°“* would not change if V were chosen to 
g g 
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be a spherical volume with its boundary S' residing in the far zone of the entire object. Then it is 
straightforward to derive that 


we _ Kya yr as Ky(n'™ or" ae Ky3(a'™)U®° a maa yr. (138) 


ws = I, dr [Z,,(F, wim’ +Zoe ne jo + Z3(t, a) £Zuan P| (139) 
1 


(see Sec. 13.4 of [34]). 

It is important to recognize that although the extinction and phase matrices are inherently 
far-field quantities, Eqs. (137)-(139) are valid for any volume enclosing the entire scattering ob- 
ject even if its boundary lies in the object’s near field. Of course, a trivial modification of Eqs. 
(137)}+(139) would not work for a volume enclosing only part of the scattering object, as illus- 
trated in Fig. 11b, since in this case one would need to know the specific near-field spatial distri- 
bution of the Poynting vector over S. Eqs. (137)-(139) would also not apply if the host medium 
were absorbing [179]. 


4.8. Foldy equations 


The general scattering formalism described in Subsections 4.1—4.3 applies equally to an ob- 
ject in the form of a single body and to a fixed multi-particle group. However, when the object is 
a group of touching and/or separated distinct components then it can sometimes be advantageous 
to modify the formalism by expressing the total scattered field as a vector superposition of the 
partial fields contributed by the individual components. Specifically, let us consider the scatter- 
ing by a fixed configuration of N = 2 distinct finite particles collectively occupying the interior 
region Vp according to Eq. (67) (see Fig. 6). It has been shown in [180,181] (see also Section 


11.3 of [170]) that the solution of the VIE everywhere in space can be expressed as follows: 
Br) =E™O)+ > [ arGer)-[ arFer)- BO), ew, (140) 
V; Vi 
i=l 
where the electric field vector E,(r) “exciting” particle 7 is given by 


N 
E(r)=E™(r)+ DEF), (141) 


J(#i) =1 


the E&“(r) are “particle—particle exciting field vectors” given by 
Ey“ (r) = | d°r’G(r,r’)- } ar’7(r',r’)-E,(@’), re V; (142) 
Vj Vj 


and T, is the ith-particle dyadic transition operator with respect to the common laboratory coor- 


dinate system. The T, satisfies the integral equation 
Tr, x’) = U(r)o(r — rv’) + vi(r)| ar’G(r,r’)-Fr’,r), rr e V;, (143) 
Vj 


where the U;(r) is the ith-particle potential function given by 


29 


0; Freer, 


k?[me(r)-l, re V; —_ 


Uj(r) = 
and m,(r) is the refractive index of particle 7 relative to that of the host medium. All position 


vectors originate at the common origin O of the laboratory coordinate system. Eqs. (140)-(143) 
form the system of integral so-called Foldy equations (FEs). They automatically incorporate all 
boundary conditions at individual-particle surfaces as well as the radiation condition at infinity 
and rigorously describe the scattered field everywhere in space. Comparison of Eqs. (84) and 
(143) reveals that, quite conveniently, T, is the dyadic transition operator of the ith particle in the 
absence of all the other particles. As such, it can be considered an individual optical identifier of 
the ith component of the group. 


4.9. Frequency-domain “multiple” scattering 


Ever since Heaviside’s Electromagnetic Theory {182], the concept of “multiple” scattering 
has been quite popular in discussions of electromagnetic scattering by multi-particle groups (see 
[183,184] and references therein). To demonstrate the actual nature of this concept in the context 
of frequency-domain electromagnetics [92], let us introduce the ith potential dyadic centered at 
the origin of the laboratory reference frame according to 


U.(r, 1’) = U,(r)d(r - rT (145) 
and introduce the following operator notation: 
BE = | r’B(r, 1’) « E(r’). (146) 
Iterating Eqs. (141) and (142), we have 
: N - : N N A ‘” * 
E,=E™+ ) GRE™+ DY > GPGtE™ 
i(#i) =1 j(#i) =1 W#j)=1 
N N ' ae 
x >» », GT,GT,GT,E" ++, (147) 
j(#i) =1 W(#j)=1 m(4l) =1 


whereas substituting Eq. (147) in Eq. (140) yields the following Neumann expansion of the total 
field: 


sy oY Y hehehe a. (148) 


i=l j(#i)=1 W#j)=1 


It is clear that the Neumann series is fundamentally based on the previously mentioned fact that 


T; for each i is an individual property of the ith particle computed as if this particle were alone 
rather than a member of the group. As a consequence, it has been rather common to characterize 
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Eq. (148) as the “order-of-scattering expansion” for the N-particle group. It should be remem- 
bered however that the FEs have a solution even when the corresponding Neumann series (148) 
does not converge. Numerical examples of possible divergence can be found in [185,186]. 


4.10. Far-field Foldy equations 


In principle, the FEs can be solved numerically in order to compute the field scattered by a 
fixed finite configuration of arbitrarily positioned particles. However, the solution becomes im- 
practicable quite rapidly with increasing N. One way to simplify the problem and make it tracta- 
ble is to consider a very sparse configuration by assuming that: 


e the N particles are widely separated so that each of them resides in the far zones of all the oth- 
er particles; and 


e the observation point is located in the far zone of any particle (but not necessarily in the far 
zone of the entire group). 


Specifically, if the incident electric field vector is given by Eq. (75) then the FEs imply that the 
total electric field vector is still given by the superposition (73), where the scattered electric field 
vector is now given by [34,181] 


N N 
EM (r) =)" g()| Af.) E™(R)) + ° AC, Ry) Ey | (149) 


i=l j(#i) =1 


Here, A, (n’,n) is the far-zone scattering dyadic of particle i centered at the particle’s own origin 
O, (Subsection 4.4); 


g(r) = i). (150) 
' 


and the vectors Ej satisfy the following system of equations: 


N 
Kj g(Rj)A;(Rj,n")- E™(R;) + g(Ry) > A(Ry, Rj) Ej, iJ =1,...,N, JFL 
#j)=l 


(151) 


The vector notation used in Eqs. (149) and (151) is explained in Fig. 12; a hat denotes a unit vec- 
tor in the respective direction. Eqs. (149)-(151) are called the far-field FEs. It is evident indeed 
that the linear algebraic system (151) is much simpler than the initial system of integral equa- 
tions (141)—(142). 

Eqs. (73) and (149) show that the total field at any observation point located sufficiently far 
from any particle in the sparse multi-particle configuration is a superposition of the incident 
plane wave and N partial spherical wavelets centered at the N particles. Importantly, the observa- 
tion point r does not have to be in the far zone of the entire N-particle group: it can be anywhere 
in space (e.g., between particles i andj in Fig. 12) as long as it is in the far zone of any particle 
entering the group. The total scattered field (149) is not, in general, a transverse electromagnetic 
wave. It becomes a transverse wave only in the far zone of the entire N-particle configuration 
defined by the criteria (89)—(91), where a is the radius of the smallest sphere encompassing all N 
particles. 
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4.11. Dyadic correlation function and Poynting—Stokes tensor 


We have already mentioned in Section 2 that it can be possible in some cases to derive an 
analytical expression for an optical observable that is explicitly devoid of the electromagnetic 
field. Sometimes this expression is a closed-form equation solving which can serve as a highly 
efficient means of calculating the optical observable directly, without the prior detailed computa- 
tion of the electromagnetic field itself. 

The general methodology enabling one to bypass an explicit use of the electromagnetic 
field is well exemplified by the far-field formulas (95), (98), and (117), in which case the ob- 
servable in question is the 4-element Stokes column vector. However, this observable can be de- 
fined only for a transverse (i.e., plane or spherical) electromagnetic wave, whereas the total elec- 
tromagnetic field in the near zone of any object (e.g., at any observation point inside a cloud of 
particles) is never a transverse wave. Furthermore, the Stokes column vector contains no explicit 
information on the direction of the Poynting vector and cannot be used in situations when this 
direction is not known a priori. 

The Poynting vector is another optical observable extensively discussed in the preceding 
sections. Its obvious analytical limitation is that different pairs of electric and magnetic field vec- 
tors can yield the same Poynting vector. This implies that the Poynting vector cannot be used to 
describe the phenomenon of electromagnetic scattering by, for example, expressing the Poynting 
vector of the scattered field in that of the incident field. In other words, the Poynting vector does 
not carry sufficient information about the participating electric and magnetic fields and, in par- 
ticular, carries no information about the polarization state of a transverse electromagnetic field. 

It is therefore highly desirable to introduce an alternative quantity that: 


e can be defined for any electromagnetic field; 
e has the dimension of electromagnetic energy flux; and 


e enables a complete and self-contained description of electromagnetic scattering in the context 
of practical optical analysis. 


It has been shown in [34,187] that a rather general quantity satisfying all these require- 
ments is the so-called dyadic correlation function involving electric field vectors at two different 
points in space: 


Cr’, r) = E(r’) ® [E(n)]’. (152) 
This quantity along with the Maxwell curl equations (68) and (69) can be used to compute other 
observables, including those involving both the electric and the magnetic field vector. An im- 
portant example is the so-called Poynting—Stokes tensor defined as the dyadic product of the 
magnetic and complex-conjugated electric field vectors taken at the same point in space: 

P(r) = + He) @ E*(r). (153) 
Indeed, it is easily verified that 

Bw) = ——1V, xen, (154) 

21MLU aa 


where the subscript r’ means that the V operator acts only on E(r’). Unlike the Stokes parame- 
ters, this quantity is applicable to any electromagnetic field (e.g., the near field of a scattering 
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object) and not just to a transverse electromagnetic wave. Furthermore, unlike the Poynting vec- 
tor, the Poynting—Stokes tensor preserves all the information about the participating electric and 
magnetic fields that gets lost upon taking the vector product of these fields. We will see in the 
following sections that owing to their generality, the dyadic correlation function and the Poyn- 
ting—Stokes tensor enable the derivation of compact closed-form analytical formulas and equa- 
tions directly describing electromagnetic scattering in terms of optical observables. 

It is straightforward to see that with respect to the Poynting—Stokes tensor, the Poynting 
vector and the Stokes parameters are derivative quantities. Indeed, we have in general 


S(r) = Re{[P5(r) — Py (W)]% + Po) — POY + [Px (0) — Py ()] 8}, (155) 
where x, y, and z are the unit vectors of a right-handed Cartesian coordinate system. In the 
case of a transverse electromagnetic wave, the Stokes column vector is given by 
Poo (t) — Pog (X) 
Pog (t) + Poo (r) 
* * 
Poo (r) - Poo (r) 
| iLPoe () + Poy (¥)] | 


I(r) = (156) 


where 6 and @ are the zenith and azimuth angles defining the local direction n of wave propa- 


gation such that n = 0 x @. 
According to Eqs. (73), (75), and (83), the total field can be expressed as 


E(r) = S-(r,n™)- ES’, re R, (157) 


inc 


where 8 e(r,n"’) is a transformation dyadic independent of Ej° and given by 


Se (r,n') = exp(ikn'™ -r)I + | d°r’Gr,r’) - | rT (r’, r”) exp(ikin'™ - r’). (158) 
Vint Vint 


It follows from Eqs. (68), (69), (74), and (76) that a relationship similar to Eq. (157) must exist 
for the magnetic field as well: 


H(r) = Sy(r,n™)- He, re R, (159) 

where the transformation dyadic $,,(r, '") is independent of Hj and is given by 

84 (r,n™) = rad x 3. (r, ni”) x ni, (160) 
l 

Then we have for the Poynting—Stokes tensor of the total field [34]: 

P(r) = S a(r, n°) 7 pire . (Ser, ny, re KR, (161) 

where 

p= 5 fii @ [ES] (162) 


is the Poynting—Stokes tensor calculated separately for the plane-wave incident field. 
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We will see below that Eq. (161) is a general template for many closed-form relationships 
between observable characteristics of the incident and total fields. Importantly, this formula 


demonstrates that the elements of the tensor P(r) generally depend on all the elements of the 


tensor P'"*. In other words, as we have already mentioned, the complex Poynting vector of the 
total field cannot be uniquely expressed in that of the incident field. This implies that the wide- 
spread characterization of electromagnetic scattering as causing the transformation of the intensi- 
ty of the incident light into that of the scattered light is fundamentally wrong. 

Eq. (161) implies the existence of a linear (but not necessarily tensorial) operator express- 
ing the Poynting—Stokes tensor of the total field in that of the incident plane-wave field [34]. We 


aS 


will denote this operator by S$ and write symbolically: 
P(r) = S(r,n™)P™, (163) 


The reader may find it instructive to rewrite Eq. (163) in the matrix form with respect to the Car- 
tesian laboratory coordinate system and thereby express the elements of the 9 x9 matrix repre- 
senting the operator S$ in terms of the elements of the 3x3 matrices representing the dyadics 
Sp and Sy. Analogously, we can write 


ps? (r) = Gsca (r, nin \P™. (164) 
where 
Poe : HS (r) ®[E(r) J" (165) 


is the Poynting—Stokes tensor calculated separately for the scattered field. 

Eqs. (163) and (164) represent a remarkably compact yet general way of describing elec- 
tromagnetic scattering in terms of optical observables rather than macroscopic field vectors. As 
such, they will be central to the following discussion, especially when it comes to the scattering 
of quasi-monochromatic fields by temporally variable objects. The reader can verify that Eqs. 
(98) and (117) are but specific coordinate-dependent manifestations of these formulas. 

Formulas analogous to Eqs. (163) and (164) can be derived for optical observables other 
than the Poynting—Stokes tensor. Each such formula serves as a linear transformer with an opti- 
cal observable of the incident electromagnetic field as the input and an optical observable of the 
total or scattered electromagnetic field as the output. Such linear transformers are essential in 
practice because of the two-layer structure of electromagnetics discussed in Section 2. 


4.12. Quasi-monochromatic scattering by a fixed object 


The formalism summarized above provides an efficient means of computing time-averaged 
macroscopic optical observables without solving explicitly the microscopic Maxwell—Lorentz 


equations. It is based, in particular, on the assumption that the complex amplitudes Ee and 


HH’ = Je, /tp a'"° x Ei® entering the solution of the standard scattering problem are inde- 
pendent of time. Let us now imagine a situation wherein these amplitudes remain constant over 
periods of time 7; such that 


T; >> Ty (averaging strategy 1) (166) 
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and 


T; > T’ (averaging strategy 2), (167) 

but fluctuate over longer time scales. In other words, 

E'(r,) = E}* (1) exp(ikn™ - r), (168) 

H(r, t) = Hi) exp(itii™ - r) = [= a’ x EP°(1) exp(ikin'™ - r), (169) 
Ho 


where significant random changes of the complex amplitude En’ (¢) occur over periods of time 


longer than 7;. The solution of the standard scattering problem for a temporal succession of 
1g (t)-values then yields a temporal succession of the total field vector values that can be 
thought of as defining time-dependent macroscopic field vectors fluctuating randomly on time 
scales longer than T;: {E(r), H(r)} => {E(r, £), Hr, t)}. 

According to the above discussion, the quasi-instantaneous values of E(r, 1) and Hr, t) 
are postulated to define optical observables averaged over time intervals of the order of 7;. For 


example, 
_ 1 ~ = . 
S(r, t) = uaa t) x[H(r, Oy} (170) 
and 
P(r, t) = = Hic, t) @ E*(r, 2). (171) 
The corresponding averages over much longer time intervals are then calculated according to 

§ Lp aS), T 

r,t) =— t S(r,t ), >T,, 172 

Geo =F), Se) (172) 

ss 1 pity. ; 
(Pony =f dP), TT, (173) 

T 4 t-T/2 


These averages are time independent provided that S(r,f) and P(r,t) are stationary random 
processes (see, e.g., [188]). 

The random macroscopic field vectors E(r, t)and Hr, t) are traditionally said to represent 
a quasi-monochromatic macroscopic electromagnetic field. In particular, Eqs. (168) and (169) 
are said to describe a quasi-monochromatic plane electromagnetic wave (or a quasi- 
monochromatic parallel beam of light). 

Despite the inequalities (166) and (167), typical fluctuations of quasi-instantaneous opti- 
cal observables still occur too rapidly to be traced by many optical instruments. It is therefore 
postulated that the intrinsic functionality of such instruments is to record the integral of an opti- 
cal observable over an extended period of time without resolving the quasi-instantaneous values 
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of this observable explicitly.’ Thus the practical usefulness of the notion of a quasi-monochro- 
matic electromagnetic field turns out to be two-fold. First, it helps combine the simplicity of the 
frequency-domain scattering formalism with a more realistic representation of the majority of 
artificial and natural sources of the electromagnetic field. Second, it allows one to account for 
inherent limitations of typical optical devices. 

The generalization of Eqs. (161), (163), and (164) to the case of quasi-monochromatic 
scattering by a fixed object is quite eee nee 


(P(r, t))) = Si (r, ni) (PO) [Sean], re R, (174) 
(P(r, t)) = Sr, n°) P™ (2), (175) 
(P(r, 2))) = 35 r, f(a). (176) 


The quasi-monochromatic versions of the main formulas of Subsections 4.6 and 4.7 are again 
coordinate-specific manifestations of Eqs. (175) and (176): 


(Signal 1(¢)) = 5 (1* ra, Ny) = 2 ZH, AYE (H)), (177) 
r 

(Signal 2(t))) = Sy (I'"(o)) — K(a'™®)(E""(2))) + “ Za’, a) (LE (2))), (178) 
(W(t) = (Wo) — (WS ))), (179) 
(W™ ())) = Ky (H™ (EO) + Ky )(O™ (0) 

+ K,3(@™)(U"* (1) + Kya (V™ @))), (180) 
(WS (t))) = I... dé [Z,,(F, nT" (2) + Z,9 (FH K(O™ (1) 

+ Zi3(F, nm )(U'™ (t))) + Zig(F, (V(t) ], (181) 

where 


| (Eby OLEdg (OT + Eq (LE Op (OT) ] 
(Er (1) = 1 [& (E00 OLEde (OF - Exp OLEop OF) | 
2 V Ho | (£09 OLE (OT - Eo LEG OT) 
| (Egy OLE OF - Eos OLEsp OT) | 
All time averages in Eqs. (174)-(182) are independent of time. 


All results of this subsection can easily be generalized to the case of a polychromatic inci- 
dent field with quasi-monochromatic components [34]. 


(182) 


7 The ultimate justification of this postulate must be based, in particular, on the explicit QED 
treatment of light—matter interactions | 110,189]. 
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5. Electromagnetic scattering by a randomly changing macroscopic object 
5.1. Dynamic and static scattering 


So far we have been discussing electromagnetic scattering by a fixed macroscopic object. 
In the case of a randomly changing macroscopic object such as a DRM, temporal changes in par- 
ticle positions and/or physical states result in significant variations in the solution of the standard 
scattering problem even if the incident field is monochromatic. The typical time interval over 
which macroscopic quadratic and bilinear forms in the field vary significantly will be denoted by 
T,. We will assume hereinafter that 7,>>7,, T,>>T’, and T,>> 7p. 

In some cases, the temporal resolution of optical measurements is finer than 7), 1.e., is suf- 
ficient to trace the random variations in macroscopic optical observables. Such measurements 
and their theoretical simulations constitute the subject of dynamic light scattering {190,191}. In 
other cases such random variations occur too rapidly to be captured by an actual optical device 
accumulating the signal over an extended period of time. This type of measurements and their 
theoretical modeling belong to the discipline of static light scattering [34]. If T is the integration 
time of an optical measurement defining its temporal resolution, then 


To 
(7 «K<T;<T<T, (183) 
for dynamic scattering and 

Ty 

r <7, <7, <T (184) 


for static scattering. In either case the practical quantification of electromagnetic scattering by a 
stochastic macroscopic object requires, strictly speaking, repeated solutions of the standard scat- 
tering problem for temporally evolving instantaneous states of the object. 

In what follows, we will mostly discuss static scattering of monochromatic and quasi- 
monochromatic electromagnetic fields. 


5.2. Monochromatic static scattering by a randomly changing macroscopic object 


To represent an actual static measurement, quadratic and bilinear forms in the field must be 
averaged over a sufficiently long period of time 7 >>7,. In the case of monochromatic scattering, 
Eqs. (170) and (171) become 


(S(r,0))) = 5 Rel ‘)x[H(r,)))), re Vexr. (185) 


(P(r, t))) = 5 (Hier, t)®@[E(r,t)]")), re Vexr, (186) 


where the macroscopic field vectors K(r, t) and H(r, t) depend on time owing to the temporal 
changes of the scattering object, while «S(r,t))) and (P(r, t))) are time independent provided 


that S(r,) and P(r,t) are stationary random processes. Now the temporal average on the right- 


hand side of Eq. (185) or (186) cannot, in general, be expressed as a product of the individual 
averages, 
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(S(r, 1))) # ; Re[(E(r, £))) x (A(r, 0))"], (187) 


(P(r, t))) # 5 (Her, t))) ®@ (E(r, t)))*, (188) 


and must be calculated explicitly. As usual, this computation is drastically simplified by assum- 
ing ergodicity of the scattering object and the resulting ergodicity of the random processes 


S(r,¢) and P(r, 7) : 
(S(r,0)) = Ge, Y))y = f dv Ser, ¥)PCH) 

= 5 Re | d¥ K(r, ¥) x [H(r, YT PY), (189) 
(Perso) = (Pe. Y)y = | av Poe) PC) 

= : | d¥Y Hr, Y) ® [E(r, Y)] PY), (190) 
where Y defines the state of the macroscopic object (rather than that of the constituent mole- 


cules) and P(Y) is a suitable time-independent probability density function. 
Egs. (161), (163), and (164) now become 


(P(r, t))) = (S(t, ns 2) + PP -[3e(r, ns 0]"*)) 


= (Sy (r,n™;P)- P™-[Sp(rn™ YY )y, re KR, (191) 
(Peyer's) Pr =e sy )yr™. (192) 
(pi (r, t))) = Es (r, an, t))) Pm = (gi (r, nize gp Pe (193) 


while the main formulas of Subsections 4.6 and 4.7 take the form 


So 


(Signal 1(4))) = => (ZG, n°; Py", (194) 
r 
(Signal 2(1))) = Soll — CK (A; Ye + 22 Za, As YY) A, (195) 
r 
(W(t) = (W™ (t))) — (WE (0), (196) 
(WO) = (Kia; Pw + (Kir (a; P)) yO" 
+ (Ki3(n'™; P))yU™ + (Kia(a'; 2) eV", (197) 


(WO) = | dE Ziu, As Pv 1 + Zia; YY) O™ 


An 


+(Z3(F, 0"; P))eU"™ + (Zi4(F, n'"; Y))eV"™]. (198) 
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Again, all time averages in Eqs. (191)—(198) are independent of time. 


5.3. Quasi-monochromatic static scattering by a randomly changing macroscopic object 


Let us now consider the situation wherein the solution of the standard scattering problem 
fluctuates in time owing to random temporal variations of both the incident field and the macro- 
scopic object. Eqs. (163) and (164) now become 


(P(r, t)) = (3(r, ns )P™ (1), (199) 
(PS (r, t))) = (354 (x, 0"; 1) P'™(0))), (200) 


where the averages are taken over a period of time much longer than both 7; and 7,. It is rea- 


sonable to assume that morphological changes of the scattering object are completely independ- 
ent of the temporal fluctuations of the externally generated incident field. More specifically, we 


assume that S(r,n'™;7) and P'™(t) as well as 38 (r,n'°; 14) and P(t) are pairs of inde- 
pendent stationary random processes, which implies that both (P(r, t))) and (ps (r,t))) are 
independent of time and are given by 


(P(r, 0)) = (Sr, 8"; 1)))(P™"())), (201) 
(P**(r, 1) = (3 (r, *; £)))((P™*(0))). (202) 
Finally, assuming ergodicity of the scattering object, we obtain 

(P(r,0)) = (S(r, 8"; Py CP"), (203) 
(P(r, t))) = (3 (r, A; Py (P™(0))). (204) 


In other words, the time averaging of the Poynting—Stokes tensor of the incident quasi- 


monochromatic beam and the ensemble averaging of the transformation operators 5 and S™ 
are completely separated. The corresponding generalization of the main formulas of Subsections 
4.6 and 4.7 reads 


(Signal 1(¢))) = ee (ZH, SF) yw (L™(0))), (205) 


(Signal 2(1))) = Su(I™*()) — (KR; P)) wv (LO) + “3 EAA Pe I"), (206) 


(0) = (HO) — (TO), (207) 
(HO) = (Kis ye CL") + (Ks Ye KO" (0) 
+ (Kis(™ ¥))v (U™ (0) + (Kia ls Pv (V(0)), (208) 


HO) = JF KZ (FH P)yv (EO) + Zins Pw KOO) 


+ (Zi3(F, a"; P))y (U'™ (t))) + (Zia(F, ms Pe (V™ (0)))]. (209) 
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It is straightforward to generalize all results of this subsection to the case of a polychro- 
matic incident field with quasi-monochromatic components [34]. 


Let us now assume that the transformation dyadics 8 gE and 8 H as well as the electric and 
magnetic field amplitudes of the quasi-monochromatic plane-wave incident field can be repre- 
sented as sums of average and fluctuating a “P?) ae 


Se(r,n'; 2) = (Se (rn; 0) + Se(r,n™; 1) = (Se(r, a oy + Se(r,n'’; 1), (210) 
Sir(r, ms 1) = (Si (r, 0"; D)) + Bir (r, "5 1) = (Sy (rn; Y))y + Sra"), (211) 
EO" (0) = (ESO) + Ear) = Ete + EPO, (212) 

HS" (0) = (5° ()) + Hos) = Hoe + Horo), (213) 
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where r € K°, the subscript “c” stands for “coherent”, and, by definition, 


(Sr(r, A"; 1))) = 0, (214) 
(Sir, A"; 0))) = 0, (215) 
(Eb! (0)) = (216) 
(HUF (0) = 0, (217) 


where 0 is a zero dyad. Then averaging Eq. (161) over a time interval much longer than both 7; 


and 7, while assuming statistical independence of the random incident field and the random 
scattering object yields 


(P(r, t))) = (S(t, ns 2) + P(t) [Se(r,n'"; )]"")) (218a) 


= 5 Sutra n°; Py - {HS ® [EGE }"} - (Se(r, as PY] )y 


— 


+ — (Sy (rf; P))y - (HOE) ® [ENF (OT) (Ser; Py 


i 


+5 Oh@i “inc, SW). {Hire ® [Ey }- [St(r, n “inc, SY)! at 


_ 5 (Si(r.n ns Y) - CHO () ® [ESF (O)')) - [Ser as YY] yy. (218b) 


This formula can alternatively be written as 
(P(r, 1) = (Snr, a; Py « (P™()) + (Se(r, ws PI )y 
+ (Si(r, ms) (P™(0)) [Ser ns PY] yy (219a) 
= (Sy, (r, 0 Y)- (P™@) [Sen PY] yy. (219b) 


We see again that averaging the Poynting—Stokes tensor of the incident quasi-monochromatic 
field over time is completely decoupled from the ensemble averaging. This implies that to solve 
the quasi-monochromatic scattering problem, one can solve the monochromatic scattering prob- 
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lem and then make the formal substitution P'® > ((P™(0))). 


6. Effective-object methodology 


Direct computer solutions of the MMEs for morphologically complex objects can be quite 
time-consuming and in many cases impracticable. As a consequence, there has been a wide- 
spread use of phenomenological so-called effective-medium rules intended to drastically simpli- 
fy the computation (see |72—78,192—199] and references therein). Implicitly, the main idea of an 
effective-object approximation (EOA) (more commonly known as an effective-medium approx- 
imation, or EMA) is to replace a morphologically complex object, either fixed or randomly vary- 
ing in time, by a much simpler “effective” object possessing essentially the same scattering 
properties. For example, one could think of replacing the Type-1 and -2 DRMs shown in Figs. 
13a,b by homogeneous scattering objects with the same overall shape defined by the surface S, 
as shown in Figs. 13c,d. 


In terms of the transformation operators $ and 9°, one can think of the following hier- 
archy of EOAs: 


e A deterministic EOA amounts to replacing a fixed morphologically complex scattering object 
by a fixed simple “effective” object such that 


S(r, n'°) = Serr, n'*), (220) 

(rn) = SH(r, a”). (221) 
e A semi-stochastic EOA amounts to replacing a stochastic morphologically complex scattering 

object by a fixed simple “effective” object such that 

(S(r, nl; Py ~ Sew(r, n’), (222) 

(S(r, A" Py = SoH (r, Al’). (223) 


e A stochastic EOA amounts to replacing a stochastic morphologically complex scattering ob- 
ject by a stochastic simple “effective” object such that 


(Sr SY ))y = (Seer rss Perr) rer (224) 
(Sr mS Y))y = (Seer ws Poor) Yon ee) 


Note that we intentionally defined the three EOAs in terms of the linear operators 3 and 


oH acting on an optical observable rather than on the macroscopic field vectors (of course these 
definitions can be generalized to include types of optical observables other than the Poynting— 
Stokes tensor). Traditionally, however, EMAs have been introduced with the purpose of replicat- 
ing the average macroscopic field vectors rather than specific optical observables [72—78,192— 
199]. In other words, a semi-stochastic EOA would normally be introduced as a recipe for re- 
placing a stochastic morphologically complex scattering object by a fixed simple “effective” ob- 
ject such that 


(S5(r, 0; P))y ~ 32 (r, 0), (226) 
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(Su (rn; P))y =~ SH (r, 0") (227) 


in Eq. (218b). This is equivalent to defining the EMA as that replicating the “coherent field” cre- 
ated by the object. This explains why an EMA recipe would typically be formulated in terms of 
replacing an actual heterogeneous object by that made of a homogeneous material with an “ef- 
fective refractive index”. 

Unfortunately, the traditional way of formulating an EMA is somewhat limited since it 
helps determine only two out of four terms on the right-hand side of Eq. (218b). The three EOAs 
defined by Eqs. (220)-(225) are more general and useful. Furthermore, they yield automatically 
the traditional field-based EMAs in cases when the last two terms on the right-hand side of Eq. 
(218b) can be neglected. In addition, they do not rely on the contrived notion of the coherent 
field. 

To the best of our knowledge, EOAs remain unproven hypotheses since none of them has 
been derived directly from the time-domain or frequency-domain MMEs under well-defined and 
reproducible conditions. In the words of Chylek et al. written in 2000 [77], EMAs 


are not approximations in a strict mathematical sense. It is not generally possible to esti- 
mate the accuracy of a given approximation by considering the magnitude of neglected 
terms with respect to those that are kept. EMAs are often based on an ad hoc assumption 
that leads to a simplified, solvable model of a real, complicated, and usually unsolvable sit- 
uation. As a result, one is able to derive a simple or only moderately complicated prescrip- 
tion (e.g., the mixing rule) of how to calculate the average optical properties of a heteroge- 
neous composite material from the known properties and amounts of its individual compo- 
nents. Because there are no specific algebraic terms neglected and because the exact solu- 
tion of the problem is usually unknown, the accuracy of such derived effective material 
constants (effective dielectric constants or effective refractive indices of material) and the 
precise conditions for their permissible use are not easy to assess. 


Fortunately, the current availability of efficient computer solvers of the MMEs and powerful 
computer clusters makes it possible to validate EOAs numerically, at least in special cases. Re- 
cent progress in this direction will be discussed in Subsection 8.2 and Section 9. 


7. Direct computer solvers of the macroscopic Maxwell equations 


According to the preceding discussion, the main objective of the discipline of electromag- 
netic scattering by particulate objects is the computation of optical observables that can be used 
to quantify the energy budget of a macroscopic volume or the results of measurements with actu- 
al optical instruments. Alternatively, this objective can be formulated as the computation of 


quantities such as the transformation dyadics $;-(r, n°) and S;,(r,n'°) entering Eq. (161); the 
transformation operators S(r, f°) and S$°*(r,n'"°) entering Eqs. (163) and (164); the ensem- 
ble-averaged transformation operators (S(r, ni; ))y and (8°"(r, n°; ¥))y entering Eqs. 
(192), (193), (203), and (204); and their various coordinate-specific representations. Far-field 


examples of the latter are the (ensemble-averaged) extinction and phase matrices. 
Whenever possible, all these quantities should be calculated by using a direct, numerically 
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exact computer solver of the MMEs° in combination with a suitable ensemble-averaging proce- 
dure. The majority of direct solvers of the MMEs belong to one of two broad categories. Differ- 
ential-equation techniques yield the scattered field by solving the differential MMEs or the vec- 
tor wave equation in the frequency or in the time domain. Classical examples of the frequency- 
domain differential-equation techniques are the Lorenz—Mie theory for a homogeneous or radial- 
ly inhomogeneous spherical particle [1,22,200—204] and the separation of variables technique for 
a homogeneous or layered spheroid [205-211]. General differential-equation techniques applica- 
ble to an object with essentially any morphology are the frequency-domain finite-difference 
method [212,213], the finite-difference time-domain method (FDTDM) [214-218], and the 
pseudo-spectral time-domain method (PSTDM) [219,220]. In both FDTDM and PSTDM, it is 
necessary to truncate the computational domain by introducing an artificial outer boundary and 
then ensure that spurious reflections from this boundary are suppressed. This is effectively ac- 
complished by using the perfectly matched layer proposed in [221,222] and its refinements. In- 
tegral equation methods are based on the volume or surface integral counterparts of the MMEs, 
the prime examples being the frequency-domain volume integral equation method and the close- 
ly related discrete-dipole approximation (DDA) [168,223—228]. There are also hybrid techniques 
or methods that can be derived using different approaches. Furthermore, there are general formu- 
lations, such as the frequency-domain 7-matrix method (TMM) [20,24,28,35,229—235], based on 
expanding relevant electric field vectors in special mathematical functions possessing desirable 
analytical properties and then using various ad hoc techniques to compute the transition matrices 
relating the resulting columns of the expansion coefficients. 

The FDTDM and the DDA are examples of general direct solvers of the MMEs that are ra- 
ther insensitive to the object’s morphology and thus can be applied to a multi-particle group as 
well as to a compact single-body object using the same basic computer program. Some ap- 
proaches, such as the TMM, can be made much more efficient by explicitly accounting for the 
specific object’s morphology, for example, its aggregate structure. Each direct numerical solver 
of the MMEs has its own advantages and drawbacks in terms of computer memory and execu- 
tion time requirements, convergence rate, accuracy, and range of applicability. For example, the 
more traditional versions of the TMM can be less flexible than the FDTDM and the DDA in 
terms of the scattering object’s morphology, but appear to be the fastest and most accurate tech- 
niques within the range of their convergence. Further information on direct computer solvers of 
the MMEs can be found in the reviews [236,237]. 

By definition, running a direct computer solver of the MMEs yields the monochromatic 
scattering properties of a fixed object. However, the angular scattering patterns typical of a fixed 
object with a size comparable to or greater than the wavelength are typically burdened by nu- 
merous sharp maxima and minima (called speckles) that must be smoothed out to yield repre- 
sentative static-scattering results [34,92,238,239]. The necessity of repeating computations for a 
large number of realizations of a random object for the purpose of ensemble averaging still rep- 
resents a great practical challenge. For example, the computation of electromagnetic scattering 
by a DRM can require averaging over an excessive number of fixed multi-particle configura- 
tions. In this respect the advantage of the superposition TMM (STMM) [20,230,240] is the ex- 


® By definition, a direct computer solver of the MMEs is called numerically exact if it can gener- 
ate numerical results with a guaranteed number of correct decimals. The number of correct dec- 
imals may vary depending on the available computer resources and practical accuracy require- 
ments. However, all reported decimals can, in principle, be validated by modifying computer 
program settings in order to accommodate a more stringent accuracy requirement. 
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tremely efficient quasi-analytical procedure which allows one to create a fixed quasi-random N- 
particle configuration and then average relevant optical observables over all possible orientations 
of this configuration with respect to the laboratory coordinate system [241—244]. This procedure 
captures in effect an infinite continuous set of random realizations of the N-particle group, elim1- 
nates completely the notorious speckle “noise”, and yields exceedingly accurate results. 

The first studies of electromagnetic scattering by random three-dimensional multi-particle 
groups based on direct computer solutions of the MMEs [245-248] exploited the frequency- 
domain multi-sphere method [249-253] which can be considered a particular case of the STMM. 
More recently, other numerical solvers of the MMEs have been used, such as the PSTDM and its 
variations [254-258], the DDA [259-274], the FDTDM [275-277], and the hybrid finite ele- 
ment—boundary integral—characteristic basis function method [278-281]. However, the STMM 
appears to have been the most frequently used technique [268,282—328]. Studies of two- 
dimensional DRMs composed of parallel infinite cylinders have been based on the multi-cylinder 
solution of the MMEs [329-331], the PSTDM [332,333] and the FDTDM [334]. 

Given the extreme complexity of direct computer calculations of electromagnetic scattering 
by a DRM, it is imperative to characterize the accuracy of the various numerical techniques and 
certify that internal (subjective) convergence of a technique (if achieved) ensures objectively 
converged results. This can be done by comparing benchmark data generated for the same scat- 
tering object with software implementations of completely independent methods. 

As an example, let us compare far-field results obtained with five totally independent com- 
puter programs based on the STMM [244], DDA [228], invariant-imbedding TMM (II-TMM) 
[234,335], FDTDM [336], and PSTDM [337] for the same compound scatterer in the form of a 
spherical particle hosting 10 identical non-overlapping spherical inclusions (Fig. 14a). The size 
parameters of the host and the inclusions are kjR = 10 and Ayr = 2.5, respectively, where R is 
the radius of the host and r is that of the inclusions. The corresponding refractive indices relative 
to that of the infinite surrounding medium are 1.33 and 1.55 + 10.003. The coordinates of the 10 
inclusions (in units of size parameter) are listed in Table 1. It is assumed that the compound par- 
ticle is illuminated by a quasi-monochromatic plane electromagnetic wave incident in the direc- 
tion of the positive z-axis, as shown in Fig. 14b. For demonstration purposes, we define the 


4x4 dimensionless scattering matrix F(O) according to 
F(0) = = (Z(0s = @, pe; gin = 0, 0° _ oP), (228) 


where @ € [0,7] is the zenith (polar) angle, g € [0, 277) is the azimuth angle, and © is the an- 


gle between the incidence and scattering directions (i.e., the scattering angle); the ensemble aver- 
age is taken over the uniform orientation distribution of the compound scatterer; and the normal- 
ization constant C,,, 1s given by 


oun _ | dn**(Z, (0°, gpg" = 0, a _ 0; VP) yy. (229) 
4n 


Note that F is independent of ~** owing to the random orientation distribution of the scattering 
object, while Cya represents the ensemble-averaged scattering cross section 
«Ws (t))) / (1'°(t))) for the case of unpolarized incident plane-wave field (cf. Eq. (209)). It is 


easily seen that the (1,1) element of the scattering matrix F(0) (often called the phase function) 
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is normalized according to 
| dO Fi,(@) sin@ = 1. (230) 
0 


The results of our computations are tabulated in Table 2 and visualized in Figs. 15—18. Ta- 
ble 2 gives the corresponding extinction, 


Cat 
=e 231 
= (231) 
and scattering, 
as 
sca — > 232 
e mR? 4 
efficiency factors, where 
Cot = (Kin; ¥))y (233) 


is the n'"°- independent extinction cross section (W'(t))) / ((1'"*(t)y) for unpolarized incident 


light (cf. Eq. (208)). Also tabulated are the absorption efficiency factor 


Qads = Oext — Orcas (234) 

the single-scattering albedo 

w= Qsea (235) 
Oext 


and the asymmetry parameter 
(cos) = 5 | : d@ F,(@) sin © cos@. (236) 


Unlike the case with the DDA, FDTDM, and PSTDM, the averaging over orientations by 
the STMM and II-TMM computer programs is performed analytically so that the accuracy of 
computations is unaffected by simulating the uniform orientation distribution of the compound 
object by a limited set of discrete orientations. This analytical procedure also made the STMM 
and II-STM computations for the randomly oriented composite object much faster. 

The DDA simulations were performed with the code ADDA 1.2 on the computer cluster of 
the supercomputing center of the Novosibirsk State University. We used the default parameters 
of the code while controlling the discretization level by the number 7, of so-called “dipoles” 
along the particle diameter. Five values ranging from 64 to 128 were considered, corresponding 
to dipole sizes from 2/20 to 2/40, where A = 27/k, is the wavelength in the infinite surround- 
ing medium. The orientation averaging was performed with a built-in adaptive procedure which 
adjusts the number of simulated orientations of the compound object to keep the relative uncer- 
tainty in Cz caused by averaging within 10~* [228]. As a consequence, the final numerical un- 
certainty is controlled mostly by n,. To further improve the accuracy, we applied the extrapola- 
tion to the zero dipole size, as described in [338]. This procedure also provides an internal error 
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estimate which, for the majority of computed values, was adequate, i.e., was within the actual 
differences from the STMM results. Fig. 15 shows both the “raw” DDA results for n, = 64 and 
128 and the “extrapolated” ones. One can see that the quantitative agreement between the latter 
and the internally converged STMM results is quite good (the corresponding phase functions 
F\,(©) typically differ by less than 1%). 

Figs. 16-18 demonstrate a similarly impressive agreement between the STMM results and 
those obtained with the II-TMM, FDTDM, and PSTDM computer programs. Note that unlike the 
STMM, the H-TMM is based on an alternative approach to calculate the object’s T matrix which 
is more general, but can make it somewhat more cumbersome to obtain the same benchmark pre- 
cision. The FDTDM and PSTDM results have been calculated with 4/40 and 4/50 spatial grid 
sizes and exhibit expected convergence towards the STMM curves. 

This quantitative comparison of completely independent direct computer solvers of the 
MMEs obviously certifies that these five techniques can be used in reliable far-field calculations 
of electromagnetic scattering by DRMs. 

Note that when running the STMM computer program, we increasingly tightened all nu- 
merical accuracy parameters until the final results converged internally to a very high accuracy. 
While this “subjective” convergence of the STMM results does not guarantee the same “objec- 
tive” convergence, it is still expected to be a good indicator of the actual accuracy of the final 
numbers. Given the virtual absence of such benchmark numerical data in the published literature, 
we tabulate the converged STMM scattering-matrix results in Appendix A. 


8. Direct computer modeling of electromagnetic scattering by Type-1 discrete 
random media 


In this section we will discuss the results of representative calculations of electromagnetic 
scattering by Type-1 DRMs based on direct computer solutions of the MMEs. In most cases we 
will use the model of a DRM in the form of a cluster of N identical small spherical particles ran- 
domly and uniformly distributed throughout an imaginary spherical volume V with a radius R, as 
shown in Fig. 19a (after [339]). 


8.1. Far-field speckle and its suppression 

Let us first consider far-field scattering of a quasi-monochromatic plane-wave field by two 
different fixed clusters of N = 80 identical spherical particles distributed throughout an imaginary 
spherical volume with a size parameter of k,R = 40. The size parameter of the constituent spher- 
ical particles is kyr = 4 and their relative refractive index is m = 1.32. The coordinates of the 


particles forming either cluster were chosen using a random number generator, but otherwise 
they are fixed. The laboratory spherical coordinate system used to describe far-field scattering by 


either cluster is shown in Fig. 14b where, as before, the unit vectors n'° and n** specify the 
directions of incidence and scattering. The zenith and azimuth angles of the incidence direction 


are assumed to be 6’ = 0° and g'* = 0°, respectively. The incident plane-wave field is as- 
sumed to be circularly polarized in the counter-clockwise sense when looking in the direction of 
propagation, which implies that (V'™°(2))) =(7"*(t))) and (O™(2))) = (U™(t))) = 0; the dou- 
ble angular brackets denote averaging over a time interval T >> 7;. 

The two panels of Fig. 20a show the corresponding time-independent far-field angular dis- 
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tributions of the intensity (/*“*(rn*“, £))) scattered in the backward hemisphere. These intensity 
distributions were calculated using the STMM computer program described in [243] and reveal 
typical random speckle patterns. Fig. 20b shows the result obtained by averaging the scattered 
intensity over the uniform orientation distribution of the multi-particle configuration used to cre- 
ate the top panel of Fig. 20a. This orientation averaging replaces averaging over a time interval T 
>>T, and is intended to simulate averaging over uniformly random positions of all 80 particles 
by taking advantage of the efficient analytical procedure afforded by the STMM, as discussed in 
Section 7. Predictably, the average intensity pattern is rotationally symmetric with respect to the 
incidence direction and is fairly featureless, the strong and narrow backscattering peak being the 
only notable exception. 

To interpret the results of these computations, we will invoke the mathematical concept of 
ordered multi-particle sequences representing the various terms on the right-hand side of Eq. 
(148). Fig. 21a shows schematically two such sequences depicted using the blue and yellow col- 
ors. To make the discussion even more physically appealing, we will assign a cumulative phase 
to each multi-particle sequence by assuming that each particle of the sequence resides in the far 
zone of the preceding particle. For example, particle 4 of the blue sequence in Fig. 21a is in the 
far zone of particle 3, particle 3 is in the far zone of particle 2, etc. In other words, we will use 
the far-field version of the Neumann expansion (148): 


E(r, t) = E™(r, 1) + EX (r, 0) 


N 
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where we imply the notation of Fig. 12 and indicate explicitly the temporal dependence of the 
macroscopic electric field vector of the incident quasi-monochromatic plane-wave field. It is 
then easily seen that the expression for the partial electric field contributed by the blue four- 
particle sequence in Fig. 21a at the observation point includes the complex exponential factor 


exp[iki (74 + R43 + R32 + Ro; +n -R,)]. Thus the corresponding cumulative phase of the blue 
four-particle sequence is 


Optue = ki(4 + Rag + Ry + Ry, +" - Rj). (238) 


The cumulative phases of other multi-particle sequences are determined analogously. For exam- 
ple, that of the yellow three-particle sequence is given by 


Oyellow = Ai (ry + R32’ + Rory + n™ - Ry). (239) 


It is important to recognize that the very concept of the cumulative phase becomes ques- 
tionable if at least one particle of a sequence is located in the near zone of the preceding particle, 
which obviously happens in densely packed DRMs (e.g., Figs. 41 and 5). We will see however 
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that qualitative interpretations of STMM results based on the notion of the cumulative phase can 
be qualitatively instructive even in the case of random particulate volumes with substantial pack- 
ing densities. 

The origin of the far-field speckles in the two panels of Fig. 20a can now be understood by 
recognizing that in the far zone of the entire cluster the partial field due to any multi-particle se- 
quence is an outgoing transverse spherical wavelet centered at the origin of the last particle of the 
sequence. Since the distance to the far-zone observation point is much greater than the radius R 
of the imaginary particulate volume JV, all such partial wavelets at the observation point propa- 
gate in essentially the same direction given by the unit vector n°“ (Fig. 21a). The four-element 


“ sca 
n 


column (I (rn, t))) in Eq. (177) at the observation point can be directly expressed in terms 


of the elements of the scattering coherency dyadic p%*(rn%") = (ES (ra, 1) ® 


[E(-n*, f)]")) according to 


gsc ; p“ (rns) . gs + on" . p (rn) rf on" 
EI Qc 3 pr? (rns) 7 Qs ao on" . pe? (rns) ‘ ™ sca 


Lo _ Qs 7 ao (rn) : on _ ov A ep (rns) ‘ gsc ? 
i(g@ 7 p (rns) . gsc _ Qs m p (rns) . o*") 


a) 


(E* (rat, 2))) = : (240) 
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where 0° and o@* are the polar-angle and azimuth-angle unit vectors of the scattering direc- 


na 


tion such that n““ =r=0x@™. According to Eg. (237), the dyadic product 
E (rn, t) @ as (rn, ¢)]° at any moment in time is the sum of an infinite number of terms, 


each describing the result of interference of two spherical wavelets centered at the end particles 
of two particle sequences. 

Fig. 21a exemplifies one such pair. If the interference of the corresponding pair of spherical 
wavelets at the observation point is constructive (destructive) then it serves to increase (decrease) 


the total intensity scattered in the direction n**. The result of the interference depends largely on 
the phase difference A = Opine — Oyellow given by 


A= ki(% + Ry + R37 + Roy + nv’ *R) -—ry — Ryy — Ry - ni’ . Ry). (241) 


The total scattered intensity in the far zone of the particulate volume is the sum of the interfer- 
ence results contributed by all possible pairs of particle sequences. The minimal angular width of 
such interference maxima and minima is proportional to 1/A,R , while their number grows rapid- 


ly with NV. These two factors explain the typical spotty appearance of the scattering patterns in 
Fig. 20a. 

It is sometimes asserted that a speckle pattern can be caused only by monochromatic inci- 
dent light, for example by that generated by a continuous laser. In actuality, however, all one 
needs in order to observe speckles is a fixed scattering object illuminated by a quasi- 
monochromatic plane-wave field. 

The two panels of Fig. 20a exemplify the variability of the quasi-instantaneous speckle pat- 
terns that can be expected of a temporally changing DRM. After the quasi-instantaneous speckle 
patterns have been computed or measured for a representative set of evolving states of a DRM, 
one can choose to 
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e analyze the statistical information content of differences between the individual speckle pat- 
terns; or 


e apply an averaging procedure, thereby isolating the static component of the speckle patterns. 


We have already mentioned that these two approaches are known as dynamic and static light 
scattering. 


8.2. Static scattering by Type-I discrete random media 


In what follows, we simulate ensemble-averaged light-scattering characteristics of an imag- 
inary spherical volume randomly and uniformly filled with identical particles by creating only 
one random JN-particle configuration and then averaging over all possible orientations of this 
configuration with respect to the laboratory coordinate system. The fidelity of this approach will 
be analyzed later in this subsection. 

We have already seen in Fig. 20 that averaging over the equiprobable orientation distribu- 
tion of an 80-particle configuration effectively eliminates the speckle pattern and yields the com- 
bination of a smooth background and a notable backscattering peak. It turns out that the exist- 
ence of both features can be explained qualitatively by using the notion of the cumulative phase 
of a multi-particle sequence introduced above. Specifically, each far-field speckle element can be 
thought of as being the result of constructive or destructive interference of two wavelets contrib- 
uted by specific multi-particle sequences, such as those shown in Fig. 21a. The phase difference 
(241) evaluated at the far-zone observation point changes randomly as the particles move, so that 
the average result of the interference is zero. However, we will demonstrate below that certain 
classes of wavelet pairs interfere constructively irrespective of particle positions and thereby are 
responsible for the residual scattering pattern. 


ic 


Let us make a simplifying assumption that g** = g'™ and define the scattering direction 


in terms of the scattering angle © = 0°“. Then scattering in the far zone can be conveniently 
described in terms of the dimensionless 4 x 4 scattering matrix (228). Numerous STMM compu- 
tations have demonstrated that the elements populating the upper right and lower left 2 x 2 
blocks of this matrix are negligibly small compared to the other elements, which is an expected 
result of averaging over the equiprobable orientation distribution of a multi-particle group cou- 
pled with sufficient uniformity of the initial particle positions throughout the scattering volume 
(cf. Table A.1). Specifically, the scattering matrix has the following typical structure: 


Fi) Fi) 0 0 
F(0) = Fy\(O) Fy (0) - 0 _ 0 , (242) 
0 F33(0) Fy(@) 
0 0 -Fy(O) Fiy(@) 


where the scattering matrix elements denoted by a zero are at least an order of magnitude smaller 
than the smallest nonzero element (in the absolute-value sense). Note that the relations 


F2(0) = Fy (0) and Fi; (0) = ah (©) are caused by the uniform orientation distribution of a 
multi-particle cluster. In all examples discussed below, the size parameter of the imaginary 
spherical volume filled with particles is fixed at k,R = 50, while the size parameter and relative 


refractive index of the particles are fixed at Ayr = 4 and m= 1.32. 
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The plot of the phase function F\(O) in Fig. 22 reveals several fundamental consequences 


of increasing the number of particles N in the volume. First of all, there is a strong and narrow 
forward-scattering enhancement owing to the systematically constructive interference of the 
wavelets singly scattered by the constituent particles in the exact forward direction. This feature 
is detailed in the top left-hand panel of Fig. 23 and, according to Fig. 21b, can be called forward- 
scattering localization of electromagnetic waves [287]. Indeed, the left-hand panel of Fig. 21b 
shows that the exact forward-scattering direction is unique in that the phases of the wavelets 
forward-scattered by all the individual particles in the DRM are precisely the same, irrespective 
of the specific instantaneous particle coordinates [1]. It is straightforward to show that if there 
were no multi-particle sequences, the constructive interference of these single-particle wavelets 


would cause an increase of the forward-scattering phase function F,(0) by a factor of N. The 
top left-hand panel of Fig. 23 shows that this increase does occur for N = 2, 5, and 20, but even- 
tually the /,,(0) value saturates. This behavior can be explained qualitatively by referring to a 


multi-particle interaction effect whereby particle 3 in the right-hand panel of Fig. 21b “shadows” 
particle 2 by attenuating the incident field exciting particle 2. 

The second remarkable consequence of increasing N is that the phase function at backscat- 
tering angles starts to develop a narrow peak with a maximum at © = 180° (see the top right- 
hand panel of Fig. 23). The qualitative explanation of this so-called weak localization of elec- 
tromagnetic waves’ (otherwise known as the coherent backscattering effect) is illustrated in Fig. 
21c. The blue and yellow outgoing wavelets are contributed by the same chain of n particles but 
sequenced in opposite order. The opposite sequencing is largely inconsequential owing to the 
reciprocity relation for the scattering dyadic (93a). Therefore, the two conjugate wavelets inter- 
fere at the observation point constructively or destructively mostly depending on the resulting 
phase difference between the blue and yellow sequences given by 


A = k(R, —R,)- (n'™ + 0), (243) 


’ n' then the av- 


If the observation direction n°“ is far from the exact backscattering direction — 


erage effect of this interference is zero owing to randomly varying positions of particles 1 and n. 
However, at exactly the backscattering direction the differential phase A vanishes identically for 
any n-particle chain, thereby causing the interference to be always constructive and create a 
backscattering intensity peak. 

The third obvious consequence of increasing the number of particles in the DRM is the 
progressively smooth and featureless profile of the phase function at scattering angles 
30° < © < 170°. This effect manifests itself as the “diffuse” intensity background in Fig. 20b 
and is mostly caused by another class of wavelet pairs illustrated in Fig. 21d. In this case the 
wavelet caused by the yellow sequence of n particles is the same as that caused by the blue se- 
quence and thus “interferes with itself.” Since the self-interference is always constructive irre- 
spective of the specific chain of particles owing to the identity A = 0, the positive contribution 
of this class of wavelet pairs survives the ensemble averaging for any incidence and scattering 
directions. The qualitative explanation of the progressive smoothness of the phase-function 
curves with increasing N in Fig. 22 is that the side-scattered intensity is averaged over the contri- 


° Note that the frequently used term “weak localization of photons” is thoroughly inappropriate 
since it refers to an interference phenomenon that is purely classical and has nothing to do with 
QED photons. 
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butions from the rapidly increasing number of multi-particle chains. 
The bottom left-hand panel of Fig. 22 shows that the most prominent effect of increasing NV 
on the ratio — F5)\(0) j F\,(®) is to smooth out the low-frequency oscillations in the single- 


sphere curve and, on average, to make this ratio more neutral. This implies that the main contri- 


sca__A inc 


bution to the second Stokes parameter of the scattered light, ((O""(7m'”, t))), comes from single- 
particle chains, whereas the contributions from many-particle chains are largely randomized. 

A fundamental property of the ratio Fy (@)/ F,,(©) is that it is identically equal to unity 
for scattering by a single sphere [1,34]. Therefore, the rapidly increasing deviation of this ratio 
from unity for N 25 in Fig. 22 can also be attributed to multi-particle chains. Similarly, 
F;;(©) = Fi4(@) and F3(180°)/ Fi,(180°) = 1 for scattering by a single spherical particle, but 
the cumulative contribution from multi-particle chains in particulate volumes with N = 5 cause 


rapidly growing violations of these identities. 
If the incident plane-wave field is polarized linearly in the xz-plane then the angular distri- 
bution of the corresponding cross-polarized scattered intensity is defined by 


1[Fi() = Fy(0)]. This quantity is plotted in Fig. 23 along with the quantity 
1[ (0) + Fi4(0)] defining the same-helicity scattered intensity for the case of the incident 


plane-wave field polarized circularly in the counterclockwise direction when looking in the di- 


rection of the unit vector n'"® Both quantities provide the most definitive demonstration of the 
onset of weak localization with increasing N. Indeed, the corresponding single-particle curves 
show no backscattering enhancement whatsoever, so the backscattering peaks that develop with 
increasing N (and thus with growing contributions from multi-particle chains) can be attributed 
unequivocally to weak localization. 

Fig. 23 also depicts the angular profiles of the linear and circular polarization ratios defined 


as 
Fi, (©) - F(@) 

O)=a = = 244 
MUO) = 6) + 2Fix() + Fal) “4 
and 
op Fi (O) + F44(®) (245) 


~ -Fii(@) — Fus(O) 


respectively. The first quantity pertains to the case of a linearly polarized plane-wave incident 
field and is the ratio of the cross-polarized and co-polarized scattered intensities. The second 
quantity is relevant to the case of a circularly polarized plane-wave incident field and is the ratio 
of the same-helicity and opposite-helicity scattered intensities [34]. Fig. 23 demonstrates that the 
contribution from multi-particle chains serves to increase significantly the background deviations 
of both polarization ratios from zero, while weak localization causes pronounced backscattering 
peaks in the yu and wc angular profiles. 

Let us now examine whether it was indeed appropriate to calculate each ensemble- 
averaged scattering pattern in Figs. 22 and 23 by averaging over orientations of only one quasi- 
random J-particle configuration. We essentially assumed that the results thus obtained would be 
statistically representative of the average over all possible realizations of the N-particle group, at 
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least for large N. The correctness of this assumption is confirmed by Fig. 24 computed for two 
different realizations of a random 200-particle group populating a 4,R = 50 imaginary spherical 
volume. The refractive index of the identical kjr = 4 particles is again 1.32. The reader can see 
that although the two sets of initial coordinates of the 200 particles were quite different, averag- 
ing over all orientations of each configuration yielded virtually indistinguishable results. 

A more subtle and less ubiquitous manifestation of coherent backscattering can be exhibit- 
ed by a DRM populated by quasi-Rayleigh particles with sizes significantly smaller than the 
wavelength [296]. Fig. 25 depicts the ratio — F>;(O) / F\,(©) for a spherical particulate volume 
with k,R = 31 populated by N = 1, ..., 1875 identical spherical particles with kyr = 2 and m= 
1.31 [243]. It can be seen that unlike the — Fy, (0)/F,() trend in Fig. 22, the increase of NV 
first to 75 and then to 750 causes the onset and swift growth of a new feature not exhibited by the 
N= 1 curve. This narrow asymmetric minimum at backscattering angles was called the polariza- 
tion opposition effect [340]. 

Like other manifestations of coherent backscattering, the polarization opposition effect is 
caused by pairs of multi-particle sequences exemplified by Fig. 21c. A qualitative interpretation 
of this specific feature is shown in Fig. 21g using simple two-particle sequences [341]. Particles 
1-4 lie in a plane normal to the incidence direction and are assumed to have sizes significantly 
smaller than the wavelength. Particles 1 and 2 lie in the scattering plane (defined again as the 
plane through the illumination and observation directions), while the line through particles 3 and 
4 is normal to this plane. If the incident quasi-monochromatic plane-wave field is unpolarized 
then both magenta sequences contribute scattered light polarized negatively with respect to the 
scattering plane (i.e., having positive values of the Stokes parameter Q), whereas both blue se- 
quences contribute positively polarized scattered light (i.e., having negative values of the Stokes 
parameter Q). The phase difference between the conjugate magenta sequences is identically 
equal to zero, while that between the blue sequences is zero when the angle a = 180° — © (tra- 
ditionally called the phase angle) is zero, but oscillates rapidly with increasing a. Therefore, on 
average, weak localization will enhance the negatively polarized scattering contributions over a 
wider range of phase angles than the positively polarized contributions. The result is the polari- 
zation opposition effect in the form of a negative polarization minimum at a small a comparable 
to the angular width of the coherent phase-function peak. 

Despite its subtlety,'° the polarization opposition effect was observed in the laboratory 
much earlier than the more ubiquitous backscattering intensity peak. Fig. 26 shows polarization 
measurements by Lyot [342] for a particulate surface obtained by burning a tape of magnesium 
under a glass plate until the deposit on the plate was completely opaque. Lyot described the ob- 
served phase curve of polarization as “puzzling” and tentatively attributed it to the very small 
size of magnesia grains. Lyot’s results were recently reproduced and supplemented by photomet- 
ric measurements [343] (see Fig. 27). The latter revealed an equally narrow backscattering inten- 
sity peak, thereby confirming that the backscattering intensity and polarization features have 
weak localization as their common cause. The polarization opposition effect with its typically 
asymmetric angular profile was not formally identified as a manifestation of weak localization 
until 1993 [340]. However, its physical origin is precisely the same as that of the so-called azi- 
muthal asymmetry of the coherent backscattering cone observed in the late 1980s [298,344,345]. 


'0 For example, we have already pointed out that the bottom left-hand panel of Fig. 22 shows no 
signs of a sharp polarization minimum at backscattering angles emerging with increasing N. 
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It appears that Oetking [346] was the first to observe weak localization in the form of a nar- 
row intensity peak centered at the exact backscattering direction. However, neither Lyot nor 
Oetking offered a correct theoretical explanation of their laboratory results. The first theoretical 
prediction of weak localization was made by Watson [347] with a reference to a private commu- 
nication from R. Ruffine. The first deliberate laboratory demonstrations of coherent backscatter- 
ing accompanied by a correct theoretical interpretation should be credited to Kuga and Ishimaru 
[348], Tsang and Ishimaru [349], Van Albada and Lagendik [350], and Wolf and Maret [351]. 
Further references can be found in [25,38,352—354]. Remarkable manifestations of weak locali- 
zation in planetary astrophysics are discussed in [31,296,306,317,339,340,355—358]. 

Qualitatively, the effect of increasing the number of particles V in a DRM can be expected 
to be twofold. On one hand, it serves to increase the number of multi-particle sequences and 
thereby enhances such corollaries of the far-field Neumann expansion (237) as the smoothness of 
the scattered intensity at side-scattering directions and the various weak localization features at 
backscattering directions. On the other hand, it eventually yields packing density values so high 
that they cause features in the scattering patterns not implied by the far-field Neumann expan- 
sion. Therefore, the above qualitative interpretation of numerically exact STMM results can be- 
come partly or completely inadequate [304,305,318]. Fig. 25 shows that this is indeed the case: 
the black solid curve reveals a high-frequency ripple reminiscent of a homogeneous spherical 
particle with a size parameter comparable to that of the entire particulate volume. The corre- 
sponding packing density of 50% is so high that the expansion (237) along with the assumptions 
of randomness and statistical uniformity of particle positions become inapplicable. 

Despite this conclusion, the direct solutions of the MMEs displayed in Figs. 22 and 23 do 
demonstrate that the classical corollaries of the low-density limit can survive — at least in a semi- 
quantitative sense — volume packing densities reaching 30%. Such values are typical of particle 
suspensions and many particulate surfaces. 

Extensive STMM results reported in [287,296,298,304] have shown that the coherent 
backscattering peaks such as those in Fig. 23 are rounded at © = 180° owing to the finite size of 
the respective DRMs. The angular widths of the backscattering peaks and of the polarization op- 
position minimum are inversely proportional to k,R and are independent of N until the effects of 
packing density start to dominate. For the same /,R, the angular widths of the backscattering 


peaks (but not their amplitudes!) are weakly dependent on the particle size parameter and refrac- 
tive index. Mixtures of spherical particles with different size parameters or different refractive 
indices also reveal all typical manifestations of weak localization, thereby further corroborating 
the universal interference nature of this phenomenon [307]. 

In [301], the conventional orientation-averaging procedure developed in the framework of 
the STMM was generalized to include the case of illumination by a finite Gaussian beam. Exten- 
sive computations demonstrated that all scattering patterns observed in the far zone of a random 
multisphere object and their evolution with decreasing width of the incident beam can still be 
interpreted in terms of forward-scattering interference, coherent backscattering, and diffuse 
background. It was shown in particular that the increasing violation of electromagnetic reciproci- 
ty with decreasing beam width suppresses and eventually eradicates all observable manifesta- 
tions of weak localization and strongly suppresses the forward-scattering interference, while do- 
ing virtually nothing to the angular profiles of intensity and polarization at intermediate scatter- 
ing angles. 

To conclude this subsection, let us discuss the applicability of the effective-object method- 
ology introduced in Section 6 to Type-1 DRMs. Specifically, we consider the result of substitut- 
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ing an imaginary spherical volume filled with a large number of identical particles (“inclusions’’) 
by a homogeneous spherical object of the same radius, as shown in Fig. 28. Obviously, this sub- 
stitution belongs to the category of semi-stochastic EOAs. In Fig. 29, the thick gray curves depict 
the orientation-averaged far-field STMM results for an imaginary 4;R = 10 spherical volume 
populated by N = 15000 identical spherical inclusions, each having a size parameter of Ayr = 0.2 
and a refractive index of m = 1.2. For comparison, the thin black curves show the Lorenz—Mie 
results for the effective-medium counterpart of this imaginary spherical volume in the form of a 
homogeneous spherical particle with k,R = 10 and meg = 1.023115. Note that this effective re- 
fractive-index value follows from the Maxwell-Garnett effective-medium rule (EMR) [77] for 
the resulting 12% volume fraction. 

It is patently obvious from Fig. 29 that despite the extremely small size parameter of the in- 
clusions and their very large number, the Maxwell-Garnett EMR fails to reproduce the far-field 
dimensionless scattering matrix of the Type-1 particulate volume at side- and backscattering an- 
gles. In fact, the results of extensive Lorenz—Mie computations for effective refractive indices 
other than 1.023115 (not shown) revealed even worse agreement with the STMM curves. The 
likely qualitative explanation of this failure is the “bumpiness” effect wherein the discrete inclu- 
sions do not reproduce sufficiently well the perfectly smooth spherical surface of the effective 
Maxwell-Garnett scatterer responsible for the large-amplitude maxima and minima in the Lo- 
renz—Mie curves. Not surprisingly, the Maxwell-Garnett EMR reproduces the STMM extinction 
cross section and asymmetry parameter much more accurately, the corresponding ratios being 


CStM™ /CMS =1.0375 and (cos @)S™™ /(cos OM = 0.9976. 


9. Direct computer modeling of static scattering by Type-2 discrete random media 


In this section, we discuss the results of representative far-field STMM calculations for 
Type-2 DRMs to analyze how well they can be replicated by the effective-object methodology 
(see also [258,284,321—324,327,328]). For the purposes of our analysis, a heterogeneous object 
is modeled as an actual spherical body randomly filled with N identical small spherical inclu- 
sions, as shown in Fig. 30a. Following the approach outlined in the preceding section, the statis- 
tical randomness and uniformity of the object’s interior is simulated in two steps. First, we use a 
random-number generator to create a fixed yet quasi-random and quasi-uniform configuration of 
the N inclusions, while making sure that the volumes of the inclusions do not cross the object’s 
boundary and do not overlap. Second, we average all far-zone optical observables over the equi- 
probable orientation distribution of the resulting heterogeneous object using the STMM code de- 
scribed in [244]. 

The STMM results shown in Fig. 31 are obtained by assuming that the size parameter of 
the spherical host is fixed at k;R = 12, while that of the inclusions takes on values Ayr = 0.3 


and 1. The respective numbers of the inclusions are VN = 12800 and 346, both implying the same 
p = 20% volume fraction. The refractive indices of the host and the inclusions are fixed at 
Myst = 1.33 and Min: = 1.55, respectively. For comparison, we also show the results of Lo- 
renz—Mie computations for a homogeneous spherical object with the size parameter (,R = 12 
and the refractive index mim = 1.372. This refractive index provides the best fit of the Lorenz— 
Mie scattering matrix to that calculated for the heterogeneous object with N = 12800 inclusions 
and, in fact, is very close to the value myc = 1.3728 predicted by the Maxwell-Garnett EMR for 
the given host and inclusion refractive indices and the inclusion volume fraction. Again, this 
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EMR is predicated on the replacement of the heterogeneous target by an equidimensional homo- 
geneous object with the same outer boundary, as exemplified by Figs. 30a,b, and belongs to the 
category of semi-stochastic EOAs. Since the Maxwell-Garnett effective refractive index is inde- 
pendent of Ayr, the thick gray curves in Fig. 31 represent the EMR substitution for both hetero- 
geneous objects. 

It is obvious that if the boundary of the host body is perfectly spherical then the Maxwell- 


Garnett EMR must reproduce the well-known Lorenz—Mie identity By (0) i F,,() = 1. There- 


fore, a deviation of the ratio Fos (0) i F,(0) for a heterogeneous spherical object from 100% is 
the most direct and unequivocal indicator of the numerical inaccuracy of the effective-medium 
methodology. Fig. 31 shows that the inclusion size parameter A,r = 0.3 yields Fy(0)/ F,(O) 
values hardly distinguishable from 100%, whereas the inclusion size parameter k,r = 1 causes 
an obvious failure of the EMR. 

Comparison of Figs. 29 and 31 reveals that the performance of the Maxwell-Garnett EMR 
is markedly better in the case of the Type-2 DRM, probably owing to the absence of the bumpi- 
ness effect. In fact, the nearly perfect agreement between the STMM curves for kyr = 0.3 and 
the Lorenz—Mie curves in Fig. 31 provides a convincing numerical validation of the effective- 
object hypothesis underlying the Maxwell-Garnett rule for Type-2 DRMs and should motivate 
efforts to derive this rule analytically from the MMEs. Still the STMM results for kjr = 1 in Fig. 
31 show that the range of applicability of the EMR in terms of the maximal permissible inclusion 
size parameter can be quite limited. This result should also be explained by the analytical deriva- 
tion. 

In general, the optical cross sections and the asymmetry parameter are known to be less 
sensitive functions of the object’s morphology than the elements of the scattering matrix. One 
can therefore expect a somewhat better accuracy of the Maxwell-Garnett prediction of the inte- 
gral radiometric characteristics than that of the angular scattering-matrix profiles even for rela- 
tively large inclusions. This is indeed the case, the eee ratios being very close to unity 


for both inclusion size parameters: Conc © =1,0066 and (cos@)S'™M (cos @)MS = 


0.9975 for kr =0.3 and ea (a cs =1.0209 and (cos @)S™™ [cos @)MS =0.9779 for 
kyr = 1. 

To further substantiate the effective-medium hypothesis, in Fig. 32 we show the results of 
T-matrix computations for a spherical host with k,R = 10 and mys = 1.33 randomly populated 
by two kinds of Ayr = 0.3 inclusions having refractive indices Ming; =1.45 and Mincz = 1.6. 


The number of each kind of inclusions is 4000. It is seen that the 7-matrix results can be repro- 
duced nearly perfectly by the Lorenz—Mie results for a homogeneous spherical object with 
k,R = 10 and mio = 1.37. Interestingly, almost the same refractive index (memr = 1.3696) fol- 
lows from the n-component effective-mixing rule [77]. The agreement eae the respective 


extinction cross sections and asymmetry parameters is also excellent: C, Cec ng = 0.9895 and 
(cos @)S™MM [cos o6)'™ =0.9943. 

Finally, in Figs. 33 and 34 we display the 7-matrix results for two cases when the refractive 
index of the host exceeds that of the inclusions. Specifically, mposs = 1.4 in Fig. 33 and 
Most = 1.6 in Fig. 34, while the inclusions are spherical voids with mi. = 1. The other pa- 
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rameters of both heterogeneous spherical objects are as follows: 4\R =10, Ayr = 0.3, and N= 
8000. It is seen that in these two cases, the diviations of the STMM curves for the scattering ma- 
trix elements other than the phase function from their best-fit Lorenz—Mie counterparts (corre- 
sponding to mim = 1.32 in Fig. 33 and mim = 1.472 in Fig. 34) are more noticeable than before, 
while the Fy,(0)/ F,() STMM curve in Fig. 34 signals significant problems with the very 
EMA methodology. Furthemore, the corresponding Maxwell-Garnett refractive indices 
(myc = 1.3123 and 1.4694, respectively) differ substantially from their best-fit Lorenz—Mie 
values. Yet Table 3 shows that the Maxwell-Garnett refractive indices yield more accurate pre- 
dictions of the extinction cross section and asymmetry parameter than the Lorenz—Mie refractive 
indices inferred by best-fitting the STMM scattering-matrix results in Figs. 33 and 34. Again, the 
still-to-be-developed analytical theory of the macroscopic effective-medium regime will need to 
explain all these numerically exact findings. 


10. First-order-scattering approximation 


Although using a numerically exact computer solver of the MMEs is the preferred way of 
quantifying electromagnetic scattering by a DRM, the applicability of this direct approach is still 
limited in terms of the number of constituent particles and the overall size of the particulate vol- 
ume relative to the wavelength. However, there are two well-defined and often-encountered 
kinds of Type-1 DRM which allow for an explicit use of the far-field Foldy equations discussed 
in Subsection 4.10. As a result, one can derive analytically rather simple expressions or equa- 
tions for key optical observables which provide for much more efficient computations by bypass- 
ing the calculation of the electromagnetic field itself. The particles forming either kind of DRM 
are sparsely and randomly distributed, but their number N must be sufficiently small for the first 
kind or tend to infinity for the second kind. In either case the far-field conditions (89)—-(91) do 
not apply to the whole DRM, which makes it necessary to first compute the ensemble-averaged 
Poynting—Stokes tensor and then use it to quantify the energy budget of the DRM and the read- 
ing of a near-field WCR [34]. 

Let us first consider the first kind of Type-1 DRM by assuming that: 


e N is sufficiently small and the average interparticle distance is sufficiently large that in the 
framework of the Foldy equations each particle can be considered as being “excited” only by 
the incident field; 

e the N-particle DRM is observed from a distance r much greater than any linear dimension L of 
the imaginary volume V circumscribing the DRM: 


r>>L; (246) 

e the observation point is allowed to be in the near zone of the entire DRM but is assumed to be 
distant enough to reside in the far zone of any of the N particles constituting the DRM; 

e all N particles are moving randomly and independently of each other throughout the imagi- 
nary volume V; 

e the physical states of the N particles change randomly and independently of each other as well 
as independently of the particle positions, where, as before, the physical state of a particle in- 
cludes all its physical characteristics except coordinates. 


These requirements are often satisfied in laboratory and in situ measurements of light scattering 
by tenuous collections of small particles such as those discussed in [176,359—366]. 
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According to the above assumptions, the second term on the right-hand side of Eq. (141) 
can be neglected in comparison with the first term. Let us choose the origin O of the laboratory 
coordinate system close to the geometrical center of the N-particle DRM and assume that the ob- 
servation point resides close enough to be in the near zone of the entire object yet sufficiently far 
to be in the far zone of any of the N constituent particles (Fig. 35). Eqs. (75), (85), (86), and 


(140) then imply that 

ms x exp(iki;) ~ aot 

ES(r) =) exp(ikii™ - R,) PA.) - ES. (247) 
i; 


i=l 

Let us now assume that the N-particle DRM is ergodic so that we can use Eq. (219b). Also, 
all particle positions R;, as well as all particle physical states €; (and thus the corresponding 
particle-centered scattering dyadics 4,(f,,n'™)) as functions of time are considered to be inde- 


pendent random processes. This implies that averaging over all the individual-particle physical 
states and over all the individual-particle coordinates can be performed independently: 


(dy = (CRE (248) 


To average over the individual particle coordinates, we assume that the corresponding coordinate 
probability density functions are given by 
iV if R,eV, 


P,(R,) = {6 if ReV for anyi=1,..., N. (249) 


This means that the individual positions of all the N particles throughout the entire volume V are 
mutually independent and statistically equiprobable. This is consistent with the assumption that 
the average particle packing density is sufficiently small. Finally, we assume that that the angular 


dependence of the individual particle-centered scattering dyadics is weak enough that at the large 
distance r from the DRM, 


A;(f;,8)~ A;(F, 8) for any i, (250) 


where r is the unit vector originating at O and pointing in the direction of the observation point 
r (Fig, 35). 

Let us first quantify the energy budget of the entire N-particle DRM. This entails surround- 
ing the volume V by an imaginary sphere S with a radius r much greater than the volume’s typ1- 
cal linear dimension L, as sketched in Fig. 35, and evaluating the integral 


(W(t) = —Re | Pr(S)) (251) 


The explicit derivation detailed in [34] requires two more assumptions. First, the size parameter 
of the volume V must be much greater than unity: 
kL> 1. (252) 


Second, the sum of the individual extinction cross sections of the N particles forming the DRM 
must be much smaller than the geometrical cross section of the volume V. The final result, for- 
mulated here for the general case of quasi-monochromatic scattering, is as follows: 


(W(t) = KW (t))) — KW (t))), (253) 
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where 


N 
CH (OY) = YUCK; Ee CL (O)) + (Kil E))s, (O™ (0) 


i=1 


+ (Ki3(n"™; €,))z, (U™ (t))) + (Kia(a'™®; Eye, (VO), (254) 


N 
(WS (a))) = | es SZ, fs Ee, (1° (O) + (Zirh fs E))e, (OO) 
i=l 


+ (Zi3(F, in"; Ee, (UO) + (Zig 5 EN (VO). = (255) 


In the above formulas, K(n'™;é;) and Z(r,n'°;&) are the particle-centered extinction and 
phase matrices of particle i, respectively. 

Let us now consider the electromagnetic response of the two distant polarimetric WCRs 
shown in Fig. 36, each having its optical axis centered at the volume element V. Both instru- 
ments are located in the near zone of the DRM yet sufficiently far from it so that each partial 
wavelet contributing to the right-hand side of Eq. (247) becomes locally flat by the time it reach- 
es a WCR. Furthermore, although the acceptance solid angle AQ of either WCR is very small, 
its distance r from the center of the DRM is large enough that the solid angle subtended by V, as 
viewed from the WCR, is smaller than AQ. As a result, either WCR captures all N partial wave- 
lets irrespective of particles’ locations within V, while WCR 2 also captures the incident plane 
wave. 

According to Subsection 4.5, WCR 1 integrates over its objective lens the time-averaged 
Stokes column vector of the superposition of the N quasi-plane wavelets propagating in essen- 
tially the same direction r;. Since WCR 1 does not capture the incident plane wavefront, it can 


be shown [34] that the quasi-monochromatic response of WCR 1 averaged over a sufficiently 
long period of time is given by 


ain.-la So ~ ~ A inc inc 
(Signal (1) = “> (ZEA; Ede (UO). (256) 
i=l 
The {objective lens, diaphragm} filter of WCR 2 passes the incident plane wave in addition to 
the N partial quasi-plane wavelets. As a consequence, the integration of the resulting Stokes col- 
umn vector over the entrance pupil of WCR 2 yields [34]: 


(Signal 2(1))) = Sal) — Do KG"; &))z, (UO) + “3 > (Za, Hs E))g, (UO). 


i=1 i=1 


(257) 


Eqs. (253)(257) represent the so-called first-order-scattering approximation for the Type-1 
DRM in the form of a small ergodic group of sparsely distributed particles. Comparison of these 
formulas with their far-field counterparts (205)—(209) shows that the reading of a near-zone yet 
sufficiently distant WCR can be quantified by summing up the corresponding single-particle far- 
field readings. 

A fundamental consequence of the additivity of the extinction and phase matrices in Eqs. 
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(254)-(257) is that the actual N-particle DRM is optically indistinguishable from that consisting 
of N statistically identical particles, each having the same average extinction and phase matrices 
given by 


N 
K(al; 8) = — KE). 


a (258) 
(ZA; Eye = LY ea Taare) 


ao 


The matrices (K(n'; €))z and (Z(f,n'"; €))z can be thought of as being averaged over a syn- 
thetic distribution of physical states of one particle p.(¢) derived from the N individual-particle 
distributions p, (¢,). Then Eqs. (254)-(257) take the following simplified form: 


(0) = NECK a"; E))e (I (0) + (Kiz(8*; &)) (0 (0) + (Kis; 6) KU" (O)) 
+ (Kis; eV" ())) (259) 


(WO) = NI dé (Zula DAUM OD + ali a"; E))< (O™(0)) 
+ (ZiaF As) MUO) + ZuF ass VO), (260) 


(Signal 1(¢))) = “s N(ZG, n; Ee" ()), (261) 


(Signal 2())) = Sait" (0) — NK(a"™; S)reMU(O)) +S at NZ, nT; Eye" (fy). (262) 


The principal advantage of the be iceesieers approximation is that it obviates the 
need to explicitly solve the MMEs for a statistically representative set of sparse N-particle con- 
figurations and replaces this complicated task by the much simpler task of finding the far-field 
solution of the MMEs for one isolated particle followed by averaging this solution over a repre- 
sentative distribution of particle physical states. Furthermore, there is no need to satisfy the most 
challenging requirement of the far-field approximation, viz., the inequality (91), by applying it to 
the entire volume V. 

The analytical derivation of the first-order-scattering approximation does not involve an 
explicit requirement that the N constituent particles be in the far-zones of each other. Instead, the 
most important explicit requirement leading to Eqs. (259)—-(262) is that the second term on the 
right-hand side of Eq. (141) be much smaller than the first term. However, this requirement does 
imply that the average separation between the particles must be appropriately large and their total 
number N must be sufficiently small. These qualitative criteria were analyzed using numerically 
exact STMM results in [288]. Further insight can be gained from recalling that Eq. (138) is valid 
in the near zone as well as in the far zone of a DRM. Therefore, far-field STMM computations 
based on this formula should be a good test of the accuracy of Eq. (259). Table 4 shows the val- 


ues of the ratio (W°*"(t))) FOS“ / (W**"(t)))S'™™ for an imaginary spherical volume with a size 
parameter k;R = 50 randomly filled with N identical spherical particles having a size parameter 
of Air =4 and a refractive index of m = 1.32. The incident field is assumed to be quasi- 
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monochromatic and unpolarized. Also shown are the corresponding values of the packing densi- 
ty p = N(r/ R)°. It is obvious that only packing densities of one percent or less can ensure high 
numerical accuracy of the first-order-scattering approximation. 

The main difference between the far-zone formula (98) and the near-zone formula (261) is 
that the latter completely ignores the forward-scattering interference explained in Fig. 21b and 
discussed in Subsection 8.2. Yet at side- and back-scattering angles both formulas should give 
similar results provided that the main requirements of the first-order-scattering approximation 
are met. In particular, the ratios of the elements of the phase matrix must become N-independent. 
Fig. 22 shows that this is the case only when AN is smaller than 20. According to Table 4, this 
again implies that the packing density must be less than one percent. 


11. Radiative transfer and coherent backscattering 
11.1. Radiative transfer theory 


Another analytical approach directly derivable from the MMEs is what is traditionally 
called the radiative transfer theory. In this case it is assumed that: 


e the N particles forming the Type-1 DRM (Fig. 3a) are separated widely enough that each of 
them is located in the far zones of all the other particles; 


e the observation point is located in the far zone of any particle in the group (but, in general, in 
the near zone of the entire group); 


e Nis very large: N >. 


The first assumption implies the applicability of the algebraic far-field FEs (149) and (151). Ac- 
cording to the second assumption, the total field at any observation point located sufficiently far 
from any particle in the sparse DRM is the superposition of the incident plane wave and N partial 
spherical wavelets contributed by the N particles. The observation point does not have to be in 
the far zone of the entire group and can be anywhere in space, including inside the DRM, as long 
as it resides in the far zones of all the N particles constituting the DRM (see Subsection 4.10). 

The third assumption implies that we can replace the full far-field Neumann expansion 
(237) by the much simpler so-called Twersky expansion. Indeed, the terms with 7 =i and / = j/ 
in the triple summation on the right-hand side of Eq. (237) are excluded, but the terms with / =i 
are retained. Therefore, we can decompose this summation as follows: 


> ye YY YD Y Var, (263) 


gel fel i=l j=l /=1 i=1 j=l l= 
Gai es J#il#i j#i 


where 6;; is the Kronecker delta. Higher-order summations in Eq. (237) can be decomposed 


similarly. The first group of terms on the right-hand side of Eq. (263) is contributed by “self- 
avoiding” sequences of particles, whereas the second group includes contributions from se- 
quences that involve a particle more than once. The approximation introduced by Twersky [367] 
helps simplify Eq. (237) by retaining only the terms contributed by all self-avoiding multi- 
particle sequences. In the limit N — oo the Twersky approximation accounts for the overwhelm- 
ing majority of multi-particle sequences and thus can be expected to yield asymptotically accu- 
rate results. 
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Since we are dealing with a near-field problem, the solution must be based on the calcula- 
tion of the time-averaged Poynting—Stokes tensor or, more generally, the time-averaged dyadic 
correlation function (152). Using the Twersky approximation of the Neumann expansion (237), 
the Twersky approximation for the dyadic correlation function can be formulated diagrammati- 
cally according to Fig. 37. The different terms entering the expanded expression inside the angu- 
lar brackets on the right-hand side of this equation can be classified using the notation introduced 
in Fig. 38a. In this particular case, the upper and lower multi-particle sequences involve different 
particles. However, the two multi-particle sequences can involve one or more common particles, 
as indicated in Figs. 38c—f by the dashed connectors. Moreover, if the number of common parti- 
cles in a diagram is two or more then they can enter the upper and lower sequences in the same 
order, as in Fig. 38d, or in the reverse order, as in Fig. 38e. The diagrams without crossing con- 
nectors are called ladder diagrams. Two such diagrams are exemplified by Figs. 21d,e. Fig. 38f 
gives an example of a mixed diagram wherein two common particles appear in the same order 
while two other common particles appear in the reverse order. By the very nature of the Twersky 
approximation, no particle can appear in either the upper or the lower sequence more than once. 

According to the preceding discussion, the assumption of full ergodicity of the DRM al- 


lows us to replace the calculation of the time average (C(r, ¥: t))) by the calculation of the en- 

semble average (C (rr; Y%))yp = (C (r’,r;R, €))r,é, where R denotes the complete set of particle 

coordinates and € denotes the complete set of particle physical states. This problem is still very 

complex in general, but becomes more manageable if we further assume that: 

e The position and physical state of each particle are statistically independent of each other and 
of those of all the other particles. 

e The physical states of all the particles have the same statistical characteristics. 


e The spatial distribution of the particles throughout the medium is completely random and sta- 
tistically uniform. 


e All diagrams with crossing connectors in the diagrammatic expansion of the dyadic correla- 
tion function can be ignored. This is the gist of the ladder approximation [368]. 


The subsequent analytical derivation is detailed in [34] (see also [25,187]) and is not 
dwelled upon in this Report since it contains no new concepts and is a straightforward mathemat- 
ical exercise. An important intermediate step is the emergence of the following matrix integro- 
differential equation: 


AAP 


§- VIG, @) = —mo(KGS)eT,4) + mo] dq ZG. 458) <1(r. 4) (264) 


traditionally called the radiative transfer equation (RTE). Here, m9 = N/V is the average number 
of particles per unit volume; (K(q; &)) « is the single-particle extinction matrix averaged over the 
physical states of all the N particles; (Z(q, q’; €)) « is the single-particle phase matrix, also aver- 
aged over the physical states of all the N particles constituting the DRM; and 
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T(r, q) 
ie, @ =| 2°” (265) 

U(r, q) 

V(r, q) 
is the real-valued so-called specific intensity column vector. The RTE is supplemented by the 
boundary condition 


T(r, qu)|re s = O(n — qe )(I™ (a), (266) 


where S is the boundary of the Type-1 DRM (Fig. 3a), q— is any unit vector directed into the 
volume V, and 6(s) is the solid-angle delta function. Note that Eqs. (264)-(266) are valid in the 
general case of the quasi-monochromatic plane-wave incident field (168)-(169). 

It is convenient to decompose the total specific intensity column vector into so-called co- 
herent (subscript “c”’) and diffuse (subscript “d’’) components: 


T(r, @) = 6(4" — @)I.(r) + Ta(r,q). (267) 
It is easily seen that these quantities are solutions of the following boundary-value problems: 

n’ - VI.(r) = —n(K(a"™; €)) 21. (r), (268) 
T(r) res = (I (O)), (269) 


re 2 


G- VT s(r.@) = - (KG Dear) + mof  4q(ZG, 4; De Tale, 


+ m(Z(q,n""; €))<1-(r), (270) 
Ta(r, Ga) |res = 9, (271) 


where Sj the illuminated part of the boundary S and 0 is a zero four-component column. The 
obvious solution of Eq. (268) is the straightforward matrix generalization of the famous Bouguer 
exponential attenuation law [93,94,100]: 


I.(r) = exp[—nos(K(a'"’; €)) 2 (I ()), (272) 


where s is the distance between the observation point r and Sj along the straight line parallel to 
n=. 

The solution of the RTE can be directly used to compute relevant near-field optical observ- 
ables. For example, the energy-budget problem is solved by using the following formula for the 


time-averaged local Poynting vector: 
(Str,1)) = | daa7 Cr, @). 273) 


where [ (r, q), traditionally called the specific intensity, is the first element of the specific inten- 
sity column vector (265). The reading of a polarization-sensitive WCR centered around the 
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“propagation direction” q per unit time is given by 


Signal t nT Ar Sol (1) + So AQT a (r, 0'"°) if q _ ni’, 
(Bignal 450) = Sul. daTe. 4) - a re: (274) 
ADG SpAQL g(r, q) if q# ni’, 


where it is assumed that the WCR is placed inside the DRM (Fig. 39) and, as before, AQ is the 
WCR’s acceptance solid angle. 

The implications of the derivation of Eqs. (264)-(274) directly from the MMEs are quite 
profound and are discussed in [34], while the genesis of these formulas is traced in [100]. There 
are several efficient computer solvers of the RTE [25,57,59,60,64,66,369-372| which make it 
much easier to deal with the RTE than with the MMEs. The fact that the reading of the WCR can 
be modeled theoretically by solving the RTE often makes the {WCR, RTE} combination a use- 
ful optical-characterization tool. Moreover, comparison of Eqs. (273) and (274) shows that a 
WCR can be used to measure the local time-averaged Poynting vector by integrating its signal 
over the entire range qe 4 and thereby solve the energy-budget problem experimentally. 
Needless to say, to enable such optical-characterization and energy-budget applications based on 
the radiative transfer theory, the DRM must possess the specific macro- and microphysical prop- 
erties discussed in the beginning of this subsection. 


11.2. The Tyndall effect 


It is easily seen that in the absence of the integral term on the right-hand side of Eq. (264), 
the solution I(r, q) of the RTE subject to the boundary condition (266) would reduce to 


I(r, q) = 5(a — @)I.(r). (275) 


This is equivalent to using the coherent-field approximation, i.e., to keeping only the first term 
on the right-hand side of Eq. (218b). As a consequence, the reading of the WCR in Fig. 39 would 
be nonzero only if the inward optical axis of the instrument was perfectly aligned with the inci- 
dence direction (cf. Eq. (274)). The fact that a WCR immersed in or looking at a turbid medium 
and having its axis not aligned with the incidence direction can generate a nonzero signal is ex- 
plained by the presence of the integral term in the RTE causing a non-zero diffuse specific inten- 
sity column vector Ta(r,q) and the resulting inadequacy of the coherent-field approximation. 
This optical phenomenon was first identified by John Tyndall [373,374] and is often called the 
Tyndall effect. Its physical origin can be traced all the way back to the inequalities (187) and 
(188). 

Typical manifestations of the Tyndall effect primarily caused by the last term on the right- 
hand side of Eq. (270) are shown in Figs. 30c-—e. In Fig. 30c the laser beam is “invisible” when it 
passes through the glass containing pure water but becomes “visible” (1.e., causes a nonzero re- 
action of the photographic camera) when it passes through a colloidal suspension. Similarly, the 
“solar rays” become “visible” upon scattering by haze or fog particles in Figs. 30d,e. 


11.3. Weak localization 


We have seen that a major approximation in deriving Eqs. (264)-(274) was keeping only 
the ladder component of the dyadic correlation function. An improvement could be the computa- 
tion of the so-called “cyclical” component caused by pairs of multi-particle sequences exempli- 
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fied by Figs. 21c,f. Indeed, let us again consider the scattering by a Type-1 DRM as shown 
schematically in Fig. 40. The DRM is illuminated by a quasi-monochromatic plane-wave field. It 
is straightforward to show that upon statistical averaging, the contribution to the total Poynting— 
Stokes tensor of all the diagrams of the type illustrated in Fig. 41 must vanish at near-field ob- 
servation points located either inside (observation point 1) or outside (observation point 2) the 
object. However, as discussed in Subsection 8.2, there is an exception corresponding to the situa- 
tion when the observation point is in the far zone of the entire DRM and is located within its 
“back-shadow” (observation point 3 in Fig. 40). Then the class of diagrams illustrated by Figs. 
21c,f and 41c—e makes a nonzero contribution that causes the coherent backscattering effect. 
These diagrams are called maximally crossed or cyclical [375] because they can be drawn in 
such a way that all connectors cross at one point. 

The inclusion of the cyclical diagrams makes the computation of the total Poynting—Stokes 
tensor much more involved [25] and limits the range of problems that can be solved analytically 
[305,318]. A fully analytical solution has so far been derived only for a semi-infinite layer com- 
posed of nonabsorbing Rayleigh scatterers [376]. In general, no closed-form analytical equation 
similar to the RTE has been derived for the computation of the coherent component of the total 
Poynting—Stokes tensor. As a consequence, this cyclical component is often computed using the 
direct Monte Carlo summation of the cyclical diagrams [339,377]. 


11.4. Validation of the analytical theory of radiative transfer and weak localization 


By virtue of being a direct corollary of the MMEs, the radiative transfer—-weak localization 
(RT-WL) theory contains no adjustable parameters inherent in semi-empirical and phenomeno- 
logical approaches. As such, it can unambiguously be compared with computer solutions of the 
MMEs and results of controlled laboratory experiments. This is very important, since some of 
the assumptions made earlier in this section are semi-qualitative and thus need to be clarified 
quantitatively. Indeed, the RT—WL theory is fundamentally based on the asymptotic require- 
ments op <1 and N > 1, where, as before, p is the particle packing density. The first inequality 
ensures that particle positions inside the volume are random, mutually independent, and statisti- 
cally uniform. Furthermore, in the case of particles with sizes comparable to and greater than the 
wavelength, it ensures that each particle is located in the far zones of all the other particles con- 
stituting the DRM. The second inequality allows one to ignore non-self-avoiding diagrams in the 
far-field Neumann expansion (237). The combination of these inequalities implies that the over- 
all size parameter of the DRM must be much greater than unity. While these inequalities are es- 
sential in the derivation of the RT—-WL theory from the MMEs, the derivation in and of itself 
does not yield specific numerical estimates of the largest allowable packing density and the 
smallest allowable number of particles. Such estimates can only be derived from quantitative 
comparisons of the approximate RT—WL results with numerical data obtained by either directly 
solving the MMEs or performing a detailed optical experiment on a fully characterized DRM. 

An important consequence of the analytical derivation summarized in Subsection 11.1 is 
that although Eq. (274) has been obtained while assuming that the observation point r is located 
in the near zone of the DRM, the entire volume starts to behave like a single far-field scatterer as 
r — [34]. This makes it possible to validate the RT theory (alone and in combination with the 
WL theory) using far-field STMM computations and the Monte Carlo computer simulator de- 
scribed in [339,377]. Some results of this validation [308] (see also [295,310,311]) are shown in 
Fig. 42. The computations were carried out for two models of a spherical Type-1 DRM with a 
size parameter of k,R = 40. All constituent spherical particles are identical and have the refrac- 
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tive index m = 1.31 and the size parameter kjr = 2. The number of particles and the correspond- 
ing packing density are N = 250, p = 3.125% in the left-hand column and N = 500, p = 6.25% 


in the right-hand column, where, as before, 0 = Nr? / R°. Fig. 42 displays separately the RT- 


only and the combined RT—WL results. 
The comparison in Fig. 42 leads to the following instructive conclusions: 


e Although the DRMs studied contain modest numbers of particles, the packing density deviates 
from zero significantly, and the size parameter of the DRMs is moderate, the quantitative 
agreement between the exact STMM and approximate RT—WL results is quite evident. Over- 
all, this comparison confirms the mesoscopic rooting of the RT—WL theory in the MMEs 
traced in Subsections 8.2 and 11.1. 


e A scattering-angle range where the STMM and RT results disagree fundamentally is that cor- 
responding to forward-scattering directions. This result can be explained by different ways of 
treating the effect of forward-scattering interference. Indeed, in the framework of the express- 
ly near-field RT theory, this effect is incorporated mathematically in the computation of the 
exponential attenuation rate inside the particulate volume [25,34], whereas, in the framework 
of far-field STMM computations it causes the strong and narrow interference peak discussed 
in Subsection 8.2."! 


e Outside a relatively narrow range of backscattering angles, the RT-only and the full RT-WL 
results are very close. This is consistent with the physical interpretation of weak localization 
as a backscattering interference phenomenon. 


e The RT-only results do not reproduce the backscattering peaks in the phase function F,(0) 
and in the linear and circular polarization ratios defined by Eqs. (244) and (245), as well as the 
asymmetric minimum in the ratio — Fj2(@) / F\,\(©) at backscattering angles exhibited by the 


STMM results. The inclusion of the cyclical diagrams serves to reproduce these backscatter- 
ing features very closely, which is again indicative of their weak-localization nature. 


e The residual differences between the RT—WL and the STMM results at side- and backscatter- 
ing angles decrease with decreasing packing density, which is an expected result. However, 
they persist even at packing densities as small as ~3%, possibly in part because the reduction 
of p is achieved by decreasing N and thus violating more significantly the requisite inequality 
Need, 


In another recent paper [378], the RT theory was tested against the results of a controlled 
laboratory experiment. Specifically, the results of high-accuracy measurements of the Stokes re- 
flection matrix for fully-characterized suspensions of submicrometer-sized latex particles in wa- 
ter were compared with the results of a numerically exact computer solution of the RTE based on 
the so-called adding method [57,66]. The quantitative performance of the RTE was monitored by 
increasing the volume packing density of the latex particles from 2% to 10%. The results of this 
study indicate that the RTE can be applied safely to DRMs with packing densities up to ~2%. 
Radiative-transfer results for packing densities of the order of 5% should be taken with great 
caution, while the polarized bidirectional reflectivity of suspensions with larger packing densities 


'l Note that in [312] the RT exponential extinction law was reproduced by near-field STMM 
computations. 
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cannot be accurately predicted. These conclusions are generally consistent with the results of 
[308]. 


12. Fixed particulate media 


We have seen in Section 8 that the diffuse speckle-free regime naturally develops from the 
speckle regime upon averaging optical observables over changing particle positions. Further- 
more, we have seen in Subsection 11.1 that it is the averaging over random particle coordinates 
that effectively leads to the RTE (264). In the case of a fixed particulate medium such as a pow- 
der surface, a sheet of paper, or a layer of paint, the speckle regime caused by scattering of a col- 
limated monochromatic or quasi-monochromatic beam persists and is easily detectable with a 
WCR having a sufficiently fine angular resolution [44,45,379]. 

As discussed in Section 5.3 of [238], the speckle regime can get suppressed in many practi- 
cal applications owing to the use of polychromatic sources of light, uncollimated illumination, 
and/or detectors of light integrating over a wide solid angle of scattering directions. In particular, 
it is the non-detection of speckle in such applications that has led to the widespread belief that 
the RT theory or its ad hoc modifications can be used to describe electromagnetic scattering by 
fixed particulate layers. 

It is important to recognize however that the RTE has never been derived directly from the 
MMEs by averaging optical observables over a range of incidence and/or scattering directions or 
over a finite spectral range instead of averaging over varying particle positions. Therefore, the 
only way to verify quantitatively whether a fixed particulate medium can behave optically as a 
DRM is to analyze the results of direct computer solutions of the MMEs. 

Fig. 43 shows the results of STMM computations of the dimensionless scattering matrix for 
two objects. The first one is a fixed configuration of N = 200 particles with Ayr = 4 and m = 1.32 
quasi-randomly and quasi-uniformly populating an imaginary k,R = 50 spherical volume and 
yielding a 10% packing density. In this case the scattering matrix is defined according to 


Ar 


F(0) = —— Z(0 = 0, 9 = 0,0" = 0,9" = 0), (276) 


sca 


where Cyca iS given by 
Cu _ | dn**Z, (0%, op"; 9" a 0, gn = 0), (277) 
4n 


and is depicted by thin black curves. The second object is a DRM modeled by assuming a uni- 
form orientation distribution of the first object. In this case the scattering matrix is defined by Eq. 
(228) and is depicted by thick gray curves. Consistent with the discussion in Subsection 8.1, the 
sharp large-amplitude oscillations exhibited by the thin black curves represent speckles typical of 
a fixed multi-particle configuration, whereas the smooth thick gray curves are representative of a 
DRM. 

Fig. 44 is analogous to Fig. 43, but now the scattering matrix (276) of the fixed multi- 
particle configuration computed at a single wavelength is replaced by the average over a range of 
wavelengths: 

TE An sca sca inc inc 
F(0) = ——— (Z(0™" = 0,9”* =0;0™ =0,9"" =0))aa. (278) 
(Coca) Aa 
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It is assumed that (i) the incident field is a polychromatic parallel beam with quasi- 
monochromatic components, and (11) all quasi-monochromatic components have the same Stokes 
parameters (see Section 13.6 of [34]). The spectral range AA is equal to 1/10 of the central 
wavelength, which implies that k,R ranges from 47.5 to 52.5 and kyr ranges from 3.8 to 4.2. 
The numerical integration over AA was performed using a Gaussian quadrature formula with 
100 division points. 

The comparison of Figs. 43 and 44 is quite revealing. First of all, it confirms that averaging 
the scattering matrix over a finite spectral range serves as an extremely efficient suppressor of 
speckles generated by a fixed multi-particle configuration. Second of all, it demonstrates that as a 
consequence of spectral averaging the scattering properties of the fixed multi-particle configura- 
tion become very similar to those of the “morphologically-equivalent” DRM. This result [380] is 
qualitatively consistent with Eq. (241) which shows that the phase difference between two multi- 
particle sequences can be randomized not only by changing particle positions but also by varying 
the wavelength. 

Although these conclusions should be viewed as preliminary and should be corroborated by 
further research, they appear to support the conventional belief that depending on specific meas- 
urement settings (e.g., polychromatic illumination), the notion of a DRM can often be broadened 
to encompass fixed particulate media. 


13. Concluding remarks 


The overall objective of this Report was to outline the first-principles physical framework 
of the discipline of electromagnetic scattering by a (slowly varying) DRM, formulate the result- 
ing physical and mathematical problems in maximally rigorous terms, and discuss the most ro- 
bust and well-characterized ways of addressing these problems. We intentionally focused on 
numerically exact computer solutions of the MMEs as the most reliable way of obtaining pro- 
found physical insights unavailable with phenomenological and heuristic theories. We also dis- 
cussed how the first-order-scattering approximation, the radiative transfer theory, and the theory 
of weak localization of electromagnetic waves can be derived directly from the Maxwell equa- 
tions for very specific and well-defined kinds of particulate medium.'? The main advantage of 
these numerical and analytical corollaries of the MMEs is that they obviate the need to introduce 
fictitious tunable parameters and poorly defined notions such as dependent, independent, and in- 
coherent scattering; elementary volume elements; incoherent light rays; photons as particles of 
light or blobs of electromagnetic energy without phases; and collective scattering effects. The 
whole evolution of physics has been in the direction of replacing phenomenological and heuristic 
approaches with first-principles ones. A major objective of this Report was to summarize recent 
contributions to this process. 

Consistent with this objective, we stayed away from discussing phenomenological and 
semi-empirical theories of light scattering by particulate media other than the effective-medium 
approach. As explained in [91,92], facile theories such as those described in [84—90] are inher- 
ently flawed in that they are typically devoid of primordial physical parameters of a DRM in- 
volved in the solution of the Maxwell equations and instead feature numerous artificial adjusta- 
ble parameters. As a consequence, they represent little more than a conglomerate of contrived yet 


'2 A more detailed discussion of the phenomenological origin of the radiative transfer theory and 
its recent transformation into a legitimate branch of statistical electromagnetics can be found in 
[100]. 
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enticingly simple formulas intended to provide a back-of-an-envelope solution of the profoundly 
complex scattering problem. The use of freely tunable ad hoc parameters makes these models a 
flexible interpolation tool capable of fitting almost any data. The price one has to pay for this in- 
terpolation capability is that the best-fit model usually has little (if any) physical meaning. 

We hope that this Report serves as a convincing demonstration of substantial recent pro- 
gress that has made the discipline of electromagnetic scattering by a DRM a full-fledged branch 
of physical optics (or, to use a catchy term, of “disordered photonics” [381]). In particular, direct 
computer solutions of the MMEs discussed in Section 8 and straightforward analytical deriva- 
tions reviewed in Section 11 have fully confirmed the purportedly mesoscopic origin of the theo- 
ry of radiative transfer and weak localization [353,354,382—384]. Indeed, they clearly demon- 
strate how the “macroscopic” regime of this theory emerges from the “microscopic”  particle- 
level regime of Maxwell’s electromagnetics upon averaging over random realizations of a large 
sparse multi-particle group. Both theoretical and experimental studies discussed in Subsection 
11.4 (see also [385]) have revealed the inevitable breakdown of the RT-WL regime when the 
particle packing density exceeds a certain threshold. This emphasizes the importance of efficient 
computer solvers of the MMEs which have no intrinsic limitations on packing density and, in 
combination with the ever growing power of computer clusters, should eventually facilitate the 
solution of outstanding problems of unprecedented complexity. 

Still the range of scattering problems that can be solved exactly remains limited. As a con- 
sequence, approximate theories of light scattering by DRMs will still be practiced in the foresee- 
able future to handle full-scale “real-life” problems. It is therefore imperative to use advanced 
computer solvers of the MMEs as well as controlled laboratory experiments to quantify numeri- 
cal errors of approximate approaches and understand their origin. Although further research is 
still needed to better validate such popular modeling tools as the first-order-scattering approx1- 
mation, the radiative transfer equation, the theory of weak localization, and the effective-medium 
approach, significant progress has already been achieved, as discussed in Subsections 8.2 and 
11.4 and Sections 9 and 10. 

The main subject of this Report can be characterized as the direct scattering problem, i.e., 
the calculation of electromagnetic scattering by a known, well-defined system. We have not dis- 
cussed how to solve the inverse scattering problem, i.e., determine the physical characteristics of 
a particulate object by analyzing its measured scattering and absorption properties. The vastness 
of this applied discipline obviously necessitates a separate review. Similarly left out are the 
countless specific applications of electromagnetic scattering by particulate media in various 
branches of science and technology. 

In this Report we focused on isolated particulate media. Yet there is an urgent need to con- 
sider even more complex problems involving different combinations of volume and/or surface 
scattering. Good examples would be a densely packed particulate layer bounded from below by a 
plane interface and a layer of continuous fluctuating medium hosting randomly positioned dis- 
crete particles and bounded by random rough interfaces. It is safe to say that the first-principles 
treatment of such problems is still at an early stage of development [305,318,386—389]. 

Finally we note that an essential assumption made at the very outset of this Report is that 
the infinite host medium surrounding the particles is nonabsorbing. A preliminary first-principles 
analysis of the general case of an absorbing host can be found in [179,390,391]. 


Acknowledgments 


We thank an anonymous reviewer for providing helpful comments on the initial version of 


68 


this Report. We appreciate numerous insightful discussions with Yuri Barabanenkov, Matthew 
Berg, Anatoli Borovoi, Oleg Bugaenko, Petr Chylek, Helmut Domke, Joop Hovenier, Vsevolod 
Ivanov, Michael Kahnert, George Kattawar, Nikolai Khlebtsov, Nikolai Kiselev, Thomas Kulp, 
Andrew Lacis, K.-N. Liou, Pavel Litvinov, James Lock, Daniel Mackowski, M. Pinar Mengiic, 
Karri Muinonen, Antti Penttila, Thomas Reichardt, Vera Rosenbush, Yuri Shkuratov, Viktor 
Tishkovets, Cornelis van der Mee, Gorden Videen, and Edgard Yanovitskij. We thank Alexandra 
Ivanova for providing Fig. 14a, Antti Penttila for providing Figs. 19a and 30a, and Guanglin 
Tang and Jianing Zhang for help with FDTDM and PSTDM computations displayed in Figs. 17 
and 18. M.I.M., L.L., B.C., and P.Y. acknowledge continued support from the NASA Remote 
Sensing Theory Project managed by Lucia Tsaoussi, the NASA Radiation Sciences Program 
managed by Hal Maring, and the NASA ACE Project. Some of the numerical results reported in 
this paper were obtained with the “Discover” supercomputer at the NASA Center for Climate 
Simulation. J.M.D. was supported by the National Academy of Sciences of Ukraine under the 
Main Astronomical Observatory GRAPE/GPU/GRID Computing Cluster Project. M.A.Yu. was 
supported by the Russian Science Foundation grant No. 14-15-00155. 


Appendix A. Benchmark STMM results 


Owing to the equiprobable orientation distribution, the dimensionless scattering matrix 
(232) has the following symmetric structure [1,34]: 


Fi) Fi) F3(0) Fa) 
F,(0)  Fy(@)  F3(@) Fu) 
—~F3(0) -F3(0) Fs) F,()| 
Fis) F4(@)  - FO) Fus(O) 


F(O) = (A.1) 


Table A.1 is a tabulation of the 10 independent elements of the scattering matrix computed with 
the STMM program [244] for the randomly oriented compound object shown in Fig. 14a and 
specified in Section 7. Note that this table well exemplifies Eq. (246). In Table A.2, we also 
tabulate the coefficients appearing in the expansions of the numerically most significant scatter- 


ing matrix elements in Wigner d-functions di,,(@) or, equivalently, in generalized spherical 
functions P","(cos @) =i"~"den(O) [20,25,34,66,392,393]: 


max max 


Fi(®) = > at Pin(cos @) =D) ai'din(O), (A.2) 
n=0 n=0 
Fy(O) + Fi3(®) = Sa + a )P3(cos@) = S (a! + a5 )d% (0), (A.3) 
n=0 n=0 
Fn (®) - F3@) = Ss" (a3 — a3 )Ps o(cos @) = 3 (a3 — a3)d3,_5(), (A.4) 
n=0 n=0 
Fy(®) = S* at Fii(cos @) = S\ alds(0), (A.5) 
n=0 n=0 
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Nmax Amax 


F,(@) = }) B'R2(cos @) = —) Bi'd’2(®), (A.6) 
n=0 n=0 


NAmax max 


Fs4() = D_ BIF2(cos) = -) | Bidih(O). (A.7) 
n=0 n=0 

Note that 

(cos @) = a (A.8) 


The number of nonzero terms in the expansions (A.2)-(A.7) is, strictly speaking, infinite. In 
practice, however, a finite upper summation limit 7,,,, is chosen such that the corresponding 
truncated sums differ from the respective scattering matrix elements within the requisite numeri- 
cal accuracy on the entire interval © e€ [0, z] of scattering angles. All numerical accuracy pa- 
rameters in the STMM program were increasingly tightened until the numbers in Tables A.1 and 
A.2 converged to within plus/minus a few units in the last decimals given. 


Appendix B. List of acronyms 
DDA discrete-dipole approximation 


DRM discrete random medium 

EMA effective-medium approximation 
EMR effective-medium rule 

EOA effective-object approximation 
FDTDM finite-difference time-domain method 
FEs Foldy equations 


II-TMM _invariant-imbedding 7-matrix method 
MMEs macroscopic Maxwell equations 
PSTDM pseudo-spectral time-domain method 


QED quantum electrodynamics 
RT radiative transfer 
RTE radiative transfer equation 


STMM _ superposition 7-matrix method 
TMM T-matrix method 


VIE volume integral equation 
WCR well-collimated radiometer 
WL weak localization 
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Table 1 
Cartesian coordinates of 10 spherical inclusions. 


n Xn yn Zn 

1 —0.215062 6.479603 0.616824 

2 —3.756010 —2.754431 5.549602 

3 0.650697 0.515307 —0.017826 

4 —2.364920 0.805033 —4.337800 

5 5.008396 —4.096047 1.241592 

6 4.504373 —4.444519 —0.820851 

7 7.303638 0.831435 —0.230329 

8 4.725006 5.314130 1.544096 

9 —0.219794 —7.116933 —0.691158 

10 3.957806 1.528642 —4.454259 

Table 2 

Integral optical characteristics. 

Technique Oext Oe, Oavs Oo (cos @) 
STMM! 1.9104 1.8839 0.02652 0.98612 0.62184 
II-TMM 1.9105 1.8839 0.02657 0.98609 0.62186 
DDA (extrapolated) 1.9129 1.8865 0.02634 0.98623 0.62199 
FDTDM (A/50) 1.9034 1.8769 0.02650 0.98608 = 
PSTDM (A/50) 1.9129 1.8864 0.02652 0.98614 = 


'The STMM results are expected to be accurate to plus/minus one unit in the last digits shown. 


Table 3 
Extinction cross-section and asymmetry parameter ratios. 


cu™ cum (cos @)5™M™M (cos @)5™M™M 
Mho 
ci cy ce (cos @)'M (cos OMS 
1.4 1.0755 1.0235 1.0348 1.0121 
1.6 0.9061 0.9903 0.9035 0.9926 
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Table 4 
Comparison of STMM and first-order-scattering 


approximation results. 


Tq ext FOSA 
N p AO) ee 
«(Ww (t)) 
1 0.0005 1 
D 0.0010 1.0001 
5 0.0026 1.0059 
20 0.0102 1.0336 
50 0.0256 1.4451 
100 0.0512 1.6952 
200 0.1024 2.4506 
400 0.2048 4.7936 
600 0.3072 7.0776 
Table A.1 


Elements of the normalized scattering matrix calculated with the STMM computer program. 


O(deg) Fil Fy, Fi3 Fig Fo F3 Fo P33 P34 P44 

0 49.71733 0.00000 0.00000 —0.00014 49.70581 0.00322 0.00000 49.70581 0.00000 49.69506 
10 12.42281 —0.38719 0.00002 —0.00000 12.41798 0.00155 —0.00002 12.31709 1.48848 12.31480 
20 3.82084 —0.65502 0.00002 0.00037 3.81143 0.00028 0.00001 3.66047 —0.77437 3.65477 
30 2.46558 0.32314 —0.00004 —0.00033 2.45656 —0.00048 —0.00001 2.41573 0.25041 2.41197 
40 2.96876 —0.02893 -—0.00002 0.00015 2.95631 —0.00032 0.00002 2.89553 —0.57369 2.88959 
50 1.00871 0.22695 —0.00003 —0.00019 0.99662 0.00000 —0.00016 0.89948 0.28267 0.89468 
60 1.04317 —0.09226 0.00002 —0.00008 1.03129 0.00054 —0.00002 0.98187 —0.26947 0.97817 
70 0.40002 0.19768 0.00001 —0.00012 0.39089 0.00027 0.00003 0.29308 0.03631 0.29193 
80 0.53819 —0.09156 0.00001 0.00004 0.52842 0.00018 —0.00003 0.49913 —0.11133 0.49892 
90 0.24845 0.13979 —0.00013 0.00018 0.24058 —0.00015 0.00013 0.15711 —0.03949 0.15781 
100 0.16012 —0.03083 0.00000 0.00010 0.14943 —0.00039 —0.00008 0.12980 —0.01340 0.13269 
110 0.19072 0.09828 —0.00024 -0.00043 0.17921 0.00032 —0.00024 0.10879 —0.04054 0.11226 
120 0.10523 0.01397 —0.00007 0.00000 0.08762 0.00003 —0.00023 0.02627 —0.06017 0.03571 
130 0.17251 0.04556 0.00023 —0.00022 0.14978 0.00007 0.00001 0.10869 —0.00973 0.12077 
140 0.32879 0.00751 0.00023 —0.00009 0.29739 —0.00018 —0.00050 0.04533 —0.27640 0.06741 
150 0.14684 0.02090 0.00016 0.00075 0.09466 0.00010 0.00066 0.03747 —0.04800 0.07547 
160 0.76055 0.01776 —0.00121 0.00079 0.70377 —0.00073 0.00117 0.39784 —0.52574 0.44526 
170 0.40934 0.22432 0.00055 —0.00063 0.37337 —0.00021 0.00003 0.24795 —0.04579 0.27032 
180 0.38852 0.00000 0.00000 —0.00581 0.24250 0.00000 0.00000 —0.24250 0.00000 —0.09647 
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Table A.2 
Expansion coefficients calculated with the STMM computer program. 


n ay a3 a3 a4 Bi Bs 

0 1.00000 0.00000 0.00000 0.90630 0.00000 0.00000 
1 1.86553 0.00000 0.00000 1.89662 0.00000 0.00000 
2 2.29661 3.61509 3.43653. ~—-2.21923 0.09381 0.13309 
3 1.86558 2.44923 2.45022 ~—‘1.91094 0.02746 0.00610 
4 1.94148 2.23947 ~—- 2.13993 1.92623 0.04024 0.18993 
5 1.82032 1.92938 1.86526 1.79708 0.00112 ~0.05645 
6 1.96397 2.10089 2.06841 —‘1.96043 0.05146 0.21263 
7 2.17566 2.09321 2.04990 2.15775 0.12290 0.02730 
8 2.51398 2.61833 2.58468 2.50124 0.13264 0.12198 
9 2.79420 2.70146 2.69891 2.82529 0.23022 0.11972 
10 2.80854 3.03972 ~—- 2.98293 _—s- 2.79412 0.21611 ~0.04008 
11 2.83598 2.74434 2.74189 2.88242 0.29442 ~0.20541 
12 2.56817 2.85635 ~—- 2.82965 -~—s 2.61048 0.21075 0.07325 
13 2.73807 2.57889 ~—- 2.48558 ~—- 2.70521 0.26578 ~0.20495 
14 —- 2.55909 ~—s- 2.82990 ~—s-2.84695 ~—-2.63705 0.10429 0.00955 
15 3.05048 2.82332 2.72507 3.01411 0.13466 0.19354 
16 3.04930 3.32459 3.32160 3.10138 0.02963 0.03529 
i 3.36077. 3.16079 = 3.22556 ~——- 3.51389 0.27596  -0.19771 
18 2.99934 3.37988 3.27850 ~—- 3.00480 0.46426 0.79836 
19 2.08446 2.05418 2.00679 —- 2.13122 0.04756  -0.78261 
20 1.23452 1.52855. 1.26855 —1.03331 0.17120 ~0.84357 
21 0.03859 0.01828 0.00073 ~—0.02990 0.22586 0.01528 
22 0.10947 0.12174 0.10523. ~—-0..09718 0.03328 0.02163 
23 0.03300 0.03695 0.03270 ‘0.02979 0.01255 ~0.00690 
24 0.00812 0.00913 0.00803 0.00728 0.00372 0.00150 
25 0.00171 0.00192 0.00166 0.00150 0.00091 0.00026 
26 0.00032 0.00036 0.00030 0.00027 0.00019 ~0.00004 
27 0.00005 0.00006 0.00005 0.00004 0.00004 0.00000 
28 0.00001 0.00001 0.00001 ‘0.00001 0.00001 0.00000 


29 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
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Fig. 1. Examples of manmade and natural small particles. (a) Commercial glass spheres (after 
[37]). (b) Sahara desert sand (after [38]). (c) Dry sea-salt particles (after [39]). (d) A 6-mm- 
diameter falling raindrop. (e) 40-nm-diameter gold particles (after [40]). (f) Interplanetary dust 
particle U2012C11 collected by a NASA U2 aircraft. (g) Red blood cells. 


ee 


Fig. 2. (a) Natural and (b) modeled soot fractals (after [41—43]). 
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Fig. 3. Two types of discrete random medium 


throughout an imaginary volume V. (b) Type 2: 
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. (a) Type 1: particles are randomly distributed 
particles are randomly distributed throughout a 
t from that of the surrounding infinite space. 


Fig. 4. Examples of natural discrete random media. (a) Clouds of interstellar dust, arranged in 
huge patches and tentacles, appears dark when they are silhoutted against the stars in the mid- 
plane of the galaxy NGC 891. Image taken with NASA’s Hubble Space Telescope. (b) Ghostly 
glow caused by the scattering of sunlight by the interplanetary dust cloud. (c) The dusty atmos- 
phere of the comet ISON photographed on 10 April 2013 with NASA’s Hubble Space Tele- 
scope. (d) Particulate Saturn’s rings photographed from NASA’s Cassini spacecraft. (e) Jovian 
clouds photographed from NASA’s Cassini spacecraft. (f) Thin diffuse clouds in the atmosphere 
of Mars photographed from NASA’s Opportunity rover. Cirrus (g) and liquid-water (h) clouds 
in the Earth’s atmosphere. (i) Raw milk. 
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Fig. 5. Examples of natural and manmade discrete random media. (a) Cross-section of a ~2.2-um 
highly porous natural organic-matter aerosol particle (after [51]). (b) Transmission electron mi- 
crograph of a high-impact polystyrene sample cut with an oscillating diamond knife. The large 
composite particle has a diameter of ~3 um (after [52]). (c) Backscattered electron micrograph of 
the cross section of an olefin polymer blend polished using an oscillating diamond knife at room 
temperature (after |[53]). (d) Particulate surface composed of glass microspheres. (e) Electron mi- 
crograph of a paint film formed by TiO2 particles immersed in a binder. (f) Dense coating 
formed by 30-nm Y,O3 crystals. 
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Fig. 6. Standard electromagnetic scattering problem. The fixed finite scattering object consists of 
N distinct and potentially inhomogeneous components. The shaded areas collectively represent 
the interior region Vpy7, while the unshaded exterior region Vgx7 is unbounded in all directions. 
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Fig. 7. Scattering in the far zone of the object. 
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Fig. 8. Optical scheme of a well-collimated radiometer. 


Fig. 9. Examples of well-collimated radiometers. (a) 26-in refractor of the Pulkovo Observatory. 
(b) NASA’s 34-m Goldstone radio telescope. (c) NASA’s Hubble Space Telescope. (d) Human 
eye. (e) Digital photographic camera. (f) Light scattering setup built at the University of Amster- 
dam (after [176]). (g) Gershun tube (after [177]). 


1m 


i, 
ya _ ~ 
- NN 
Vd « 
, ' ‘ \ 
Scattered spherical 
wave 


Fig. 10. The response of a polarization-sensitive well-collimated radiometer depends on the line 
of sight. 
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Fig. 11. Energy budget of a finite volume enclosing (a) the entire scattering object or (b) a part of 
the object. 
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Fig. 12. Vector notation used in the far-field Foldy equations. 
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Fig. 13. Effective-medium methodology. 


Fig. 14. (a) Model compound scatterer. (b) Scattering geometry. 
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Fig. 15. Elements of the dimensionless scattering matrix computed using the STMM and DDA 
for the randomly oriented composite object shown in Fig. 14a. The n, = 64 and n,= 128 DDA 


results are shown only in the Fy / Fi, panel. 
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Fig. 16. As in Fig. 15, but for STMM vs. U-TMM results. 
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Fig. 17. As in Fig. 15, but for STMM vs. FDTDM results. 
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Fig. 19. (a) An imaginary spherical volume populated by randomly positioned spherical parti- 


cles. (b) Angular coordinates used in Fig. 20. 


a b 


Fig. 20. (a) Angular distributions of the scattered intensity for two fixed spherical particulate 
volumes. (b) As in panel (a), but averaged over random particle positions. The gray scale is indi- 
vidually adjusted in order to maximally reveal the fine structure of each scattering pattern. Fig. 
19b shows the angular coordinates used for all three panels. 
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Fig. 21. (a) Interference origin of speckle. (b) Forward-scattering interference. (c) Interference 
origin of weak localization. (d) Interference origin of the diffuse background. (e) A pair of parti- 
cle sequences contributing to the time-averaged diffuse background. (f) A pair of particle se- 
quences contributing to time-averaged weak localization. (g) Interference origin of the polariza- 
tion opposition effect. 
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Fig. 22. Elements of the dimensionless scattering matrix computed for an imaginary 4,R = 50 


spherical volume of discrete random medium uniformly populated by N = 1, 2, ..., 600 particles 
with Ayr = 4 and m= 1.32. 
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Fig. 23. Elements of the dimensionless scattering matrix and polarization ratios computed for an 
imaginary 4;R = 50 spherical volume of discrete random medium uniformly populated by N = 
1, 2, ..., 600 particles with Ayr = 4 and m = 1.32. 
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Fig. 24. Elements of the dimensionless scattering matrix for two realizations of an imaginary 
spherical volume of discrete random medium with k,R = 50, N= 200, Ayr = 4, and m = 1.32. 
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Fig. 25. Polarization opposition effects. 


0 poo ee ee ep ee 
L om 5 
L . r . . 4 
be - ry 
=> —0.4 - | 
SS L | a| 
he a) a 
Pane eo ay “i 
Sool Fee ieee : 
Tr: eat ee ett x = ES Dh 2 : a 
L | 
412 pe Pe PP 
150 155 160 165 170 175 180 


Scattering angle (deg) 


Fig. 26. Polarization measurements for a particulate surface composed of small magnesia parti- 


cles. 
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Fig. 27. Measurements of intensity and polarization of light backscattered by a particulate sur- 


face composed of small magnesia particles. 


Fig. 28. An equidimensional homogeneous spherical particle replaces the imaginary spherical 


volume filled with a large number of identical inclusions. 
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Fig. 29. Orientation-averaged elements of the dimensionless scattering matrix for an imaginary 
spherical volume of discrete random medium with 4,R = 10, N= 15000, Ayr = 0.2, and m = 1.2. 
The thin black curves show the result of using the Maxwell-Garmett approximation. 
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Fig. 30. (a,b) Heterogeneous spherical target and its effective-medium counterpart. (c—e) Mani- 


festations of the Tyndall effect. 
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Fig. 31. Elements of the dimensionless scattering matrix for randomly heterogeneous and ho- 
mogeneous spherical objects with a fixed size parameter k,R = 12 (see text). 
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Fig. 32. Elements of the dimensionless scattering matrix for randomly heterogeneous and homo- 


geneous spherical objects with a fixed size parameter k,R = 10 (see text). 
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Fig. 33. Elements of the dimensionless scattering matrix for randomly heterogeneous and homo- 
geneous spherical objects with a fixed size parameter k,R = 10 (see text). 
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Fig. 34. Elements of the dimensionless scattering matrix for randomly heterogeneous and homo- 
geneous spherical objects with a fixed size parameter k,R = 10 (see text). 
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Fig. 35. The Type-1 DRM is composed of a small number of particles sparsely populating an 
imaginary volume V and is observed from a sufficiently large distance r. 
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Fig. 36. Near-field measurements of electromagnetic scattering by a small sparse DRM. 
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Fig. 37. The Twersky approximation for the dyadic correlation function. Each arrow denotes the 
local incident field; each dot denotes the left-multiplication by the corresponding scattering dy- 
adic; and each horizontal line denotes multiplication by the corresponding g-function (150). 
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Fig. 38. Classification of various terms entering the expanded Twersky approximation for the 
dyadic correlation function. 
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Fig. 39. A WCR placed inside the DRM. The size of the WCR is exaggerated relative to that of 
the DRM for demonstration purposes. The uniform shading is intended to emphasize that the 
constituent particles move randomly throughout the volume V during the measurement. 
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Fig. 40. Electromagnetic scattering by a sparse Type-1 DRM. The size of the DRM is exaggerat- 
ed relative to its distance from observation point 3 for demonstration purposes. 


Fig. 41. Diagrams with crossing connectors. 
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Fig. 42. Scattering by a spherical Type-1 DRM with a size parameter of 4;R = 40 and packing 
densities of 9 =3.125% and 6.250%, populated with identical spherical particles with a size 
parameter of Ajr = 2 and a refractive index of m= 1.31. The solid, dotted, and thick gray curves 


depict the STMM, RT-only, and RT—WL results, respectively. The RT phase functions are shift- 
ed downward to match the RT-WL phase functions at © = 150°. 
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Fig. 43. Elements of the dimensionless scattering matrix computed for an imaginary 4,R = 50 
spherical volume populated by N = 200 particles with Air = 4 and m = 1.32. Black curves: the 


multi-particle configuration is fixed. Gray curves: the results are averaged over the uniform ori- 
entation distribution of the multi-particle configuration. 
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Fig. 44. Elements of the dimensionless scattering matrix computed for an imaginary spherical 
volume populated by N = 200 particles with m = 1.32. Black curves: the multi-particle configura- 
tion is fixed and the results are averaged over a range of wavelengths such that k,R varies from 
47.5 to 52.5 and kyr varies from 3.8 to 4.2. Gray curves: the results are averaged over the uni- 
form orientation distribution of the multi-particle configuration at a single wavelength such that 
k,R = 50 and kyr = 4. 
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